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1. Examples of Integer Programming Problems. To explain what integer 
programming problems are, it is best to start with a few illustrative examples. 

Example 1 (The Assignment Problem). A team of n people have to carry out n 
jobs, each person carrying out exactly one job. If person i is assigned to job j, a cost 
dj is incurred. Which assignment minimises the total cost? 

To formulate this problem mathematically, we introduce decision variables, con- 
straints, and an objective function: 

Decision variables: For (i,j e [l,n] := {!,... , n}), we set 




1 if person i carries out job j, 
otherwise. 



Constraints: 

Each person does exactly one job, that is, 

n 

Y,Xij = l (iG[l,n]) 



Each job is done by exactly one person, 

n 

Y,Xij=l (j€[l,n]) 

i=l 



The decision variables are binary, 

xn € B := {0, 1}, (i,je[l,n]). 



Objective function: The total cost is Ym=i Y^j=i c ij x ij- 
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Fig. 1.1. How to assign workers to vehicles? 



Thus, the assignment problem can be modelled as follows: 



n n 

m i n XX/'' <••'■'* 

i=i j=i 



.t. X{j = 1 for i = 1, . . . , '/ 



^2x tj = l for j = 1, . . . ,n, 
Xij e B for i,j = 1, . . . , n. 



More general integer programs are mathematical problems of the following form: 
Integer Programs: 



T 

max c x 



s.t. Ax < b, (componentwise) 
x > 0, (componentwise) 

x G Z", 



where A is a matrix and 6, c are vectors with rational coefficients, and where inequal- 
ities such as Ax < 6 are to be understood as componentwise inequalities. 

Binary Integer Programs: 

T 

max c x 

X 

s.t. Ax < b 

x e B" 
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Fig. 1.2. Which items should we pack into the bag? 



Mixed Integer Programs: 

max c T x + h T y 

s.t. Ax + Gy <b 

x,y>0, 
y€ IF, 

where G is a matrix and h a vector with rational coefficients. 

Example 2 (The 0-1 Knapsack Problem). A knapsack of volume b has to be 
packed with a selection of n items, where item i has volume ai and value Ci. How 
should one pack the knapsack so as to maximise the total value of items in it? 

This situation can be modelled as follows, 

n 

max CiXi 

i=l 
n 

S.t. OiXi < b, 

i=l 

x e B n , 

where the variables are defined as follows, 




1 if item i is selected, 
otherwise. 



Example 3 (The Integer Knapsack Problem) . The setup is the same as in Ex- 
ample 2, but multiple copies of each type of item are available. 



In this case the decision variables are no longer binary but arbitrary nonnegative 

3 




Fig. 1.3. How to cut this paper roll so as to minimise waste for the given size orders? 

integers, so that the model becomes 

n 

max CiXi 
i=i 

n 

S.t. CLiXi < b, 

1=1 

x > 0, 
x G Z". 



A generalisation of this problem in which more than one knapsack occur is called 
the Cutting Stock Problem. This kind of problem plays an important role in the pa- 
per industry. Another close relative is the Bin Packing Problem that occurs in the 
shipping industry. 

Example 4 (The Set Covering Problem). In a city with m neighbourhoods 
[l,m] := {1, ...,m}, n potential locations [l,n] := {!,..., n\ for fire stations have 
been identified. Let Sj C [l,m] is the set of neighbourhoods that can be served from 
location j, and let us assume that establishing a fire station at location j incurs a cost 
Cj . Where should one set up fire stations so as to minimise total set-up costs, and so 
that every neighbourhood can be served from at least one fire station. 

It is tempting to formulate this situation as a combinatorial optimisation problem, 



TC 

However, this form of the problem does not lend itself to being approached by the 
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Fig. 1.4. Where to build fire stations so that every neighbourhood is reachable within a given 
target time? 




Fig. 1.5. Graph representation of the set covering problem from Figure 1.4- 

algorithmic methods we will develop. We therefore reformulate the problem as a binary 
integer programme: 

Let the incidence matrix A = [a^] be defined by 

I otherwise. 

We introduce indicator variables Xj (j G [1,^]) as follows, 

J 1 if location j is selected, 
Xj = < 

otherwise 



The covering constraint needs to be satisfied: At least one fire station must service 
neighbourhood i for each i, 

n 

^""^ QijXj > 1, for each i G 
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Fig. 1.6. A travelling salesman problem. 



Thus, Set Covering can be modelled as the following IP, 



n 

s.t. "^^aijXj > 1, (i G [l,m]), 
i=i 

x G 1™. 



Example 5 (The Travelling Salesman Problem (TSP)). A travelling salesman 
has to visit each of n cities exactly once and then return to the starting point. For 
each pair of cities i,j G [1, n] there is a direct air link from i to j. The directed graph 
(digraph) G = (V, E) in which the vertices are the cities and the directed edges are the 
air links between them is assumed to be a complete graph. It takes Cij hours to travel 
along edge ij from city i to city j . In which order should our salesman visit the cities 
so as to minimise the total travelling time? 

Note: it may be the case that c\j ^ Cji. If c^ = Cji for all i,j G [l,7i], then we 
speak of the symmetric travelling salesman problem (STSP), and (V,E) is considered 
an undirected graph. 

Formulation as a BIP: 
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Fig. 1.7. Subtours that don't visit all nodes need to be disallowed. 



Decision variables: For all i,j £ [1, n], 



Xij — 



1 if the tour contains edge ij, 
otherwise. 



Constraints: 

The salesman leaves city i exactly once 



^2 x ij = 1 (i = l,...,n). 

r-j^i 



He arrives in city j exactly once 



To eliminate solutions with subtours, introduce cut-set constraints: 



ies j^s 



Fig. 1.8. A facility location problem. 



Thus, the TSP can be modelled as follows, 

n n 

»=i j=i 

S.t. ^ X i] = 1 (* G I 1 '™])' 

a; e B" x ". 



Example 6 (Uncapacitated Facility Location (UFL)). A set of potential depot 
locationsN = {1, . . . , n) has been identified. A set of clients M = {1, . . . , m} is known, 
each of which buys products that could be delivered from one or several of the depot 
locations. The transportation costs for satisfying all of client i 's orders from depot j 
is Cij . If depot j is in use, a fixed cost fj > is incurred. Decide which depots to 
open and what proportion of the demand of client i to satisfy from depot j so as to 
minimise the total costs. 

Modelling UFL as a mixed integer program (MIP): 
Decision variables: 

For each pair e M x N let Xij G [0, 1] be the proportion of the demand of of 
client i satisfied from depot j . 
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To incorporate the fixed costs: For each j E N, we let yj — 1 if depot j is used, 
and yj = otherwise. 

Objective function: The total cost 

c ij x ij + E] fj Vj ■ 

i£M jeN jeN 

Constraints: 

100% of client i 's demand must be satisfied 

n 

^2xij = l fori = l, 



, m. 



The are proportions, and hence, 

< < 1 Mij. 

Note that the constraints x^ < 1 are superfluous, as they are implied by Y^j=i x ij = 1- 

If any notrivial demand is being satisfied from a depot, the depot must be open, 
i.e., yj E B, with yj = 1 if 3i E M s.t. x^ > 0. Note that this condition is not in 
the form required for the purposes of integer programming. Luckily, there exists an 
equivalent formulation 1 , 

m 

^Xij<myj forj = l,...,n. 

i=l 

Thus, we find the following MIP model for UFL: 

( ^ y) -L n x„ + „ E E + E fM 

v ' i£M jeN jeN 

s.t. ^2xij = l (i E M), 

jeN 

x ^ ^ m % c? - e iV )' 

x > 0, 
y E B™. 



Example 7 (Uncapacitated Lot-Sizing (ULS)). A factory produces a single type 
of product in periods 1, . . . ,n. Production in period t incurs a fixed cost ft > for 
switching the production line on, as well as a production cost p t per unit. Earlier 
production can be kept in storage for later sale at a cost of h t per unit during period 
t. Initial storage is empty. The client demand in period t is dt- How should one plan 
the production so as to satisfy all demands and minimise the total costs? 



1 Please study this condition in detail, as this type of constraint occurs often in integer program- 
ming! It is called a big M formulation. 
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Fig. 1.9. A lot-sizing problem. 



Decision variables: 

x t is the amount produced in period t (allowed to be fractional) . 
yt = 1 if production occurs in period t, and y t — otherwise. 
s t is the amount of product held in storage during period t. 
Objective function: 

n 

}](PtXt + ftVt + h t s t ). 
t=i 

Constraints: 

No product is ever thrown away, 

s t _x + x t — d t + s t , Vt. 

Production can only occur when the machine is switched on, 

x t < My t yt, 

where M is a sufficiently large number, e.g. the total demand. Note that this is a big 
M type constraint. 

Production and storage are nonnegative, Xt,St > for all t. 

We thus arrive at the following MIP model for ULS, 

n n n 

min ^p t x t + ^ + X! htSt 
x ' v ' s t=i t=i t=i 

s.t. s t -i + x t = d t + s t (t G [1, n]), 

x t <My t (te[l,n]), 

a; t >0 {te[l,n]), 

KeB (ie[l,n]), 
s = 0, 

s t > o (ts[i,4 
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Fig. 1.10. A problem of discrete alternatives in which a linear objective function has to be 
optimised over the union of two polyhedra. This example lends itself to the geometric understanding 
of big-M constraints. 



Example 8 (Discrete Alternatives (Disjunctions)). In scheduling and other ap- 
plications, problems of the following type occur: 

T 

mm c x 

s.t. < x < u and 
either ajx < b\ (1-1) 
or (inclusive) a^x < 62, (1-2) 

where ai are vectors and bi scalars. Here "either . . . or" means that at least one of 
conditions (1.1) and (1.2) must hold. 

A big M formulation can be used to model such problems as a MIPs: 

First, introduce extra decision variables for which of (1.1), (1.2) to impose: 




1 if (1.1) is imposed, 

if (1.1) is not imposed, 

1 if (1.2) is imposed, 

if (1.2) is not imposed. 



Note that even if a condition is not imposed, it may still hold, but when it is imposed, 
it must hold. 

Second, let M > max{maxo<a;< u ajx — bi : i — 1,2}, a bound that can easily be 
computed via LP, as we shall see later. 
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The above problem of discrete alternatives can now be reformulated as follows, 

min c T x 

s.t. aJx-bi<M(l- yi ) (i = 1,2), 
yi+V2 = l, 

< x < u 
y G B 2 . 



2. Linear Programming. An important special case of an integer program- 
ming problem is one without integrality constraints, e.g., 

maxc T x 

X 

s.t. Ax < b, 
x>0. 

Such problems are called linear programming problems (or LPs). Thus, we are not 
actually referring to them as IPs, because the absense of integrality constraints renders 
the problem much easier. 

While it is interesting to study LP theory in its own right for the huge number 
of applications that can be modelled with LPs, our main interest here is the fact that 
many algorithms designed to solve general IPs rely on the concept of an LP relaxation: 

Definition 2.1. Consider the IP problem 

(IP) z* = maxc T x 

X 

s.t. Ax < b, 

x G 

By giving up on the integrality constraints Xi E Z, one obtains an LP, 

(LP) z = maxc T x 

X 

s.t. Ax < b, 

x > 0, 

called the LP relaxation of Problem (IP). 

To understand the geometric meaning of the LP relaxation, note that giving up 
on the integrality constraints has two effects on the feasible set & (the set of decision 
vectors x that satisfy the constraints of the problem), 

• & becomes larger, 

• & becomes convex. 



Proposition 2.2. The consequence of the first effect is that z > z* . 
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Fig. 2.1. The polyhedron is the enlarged ("relaxed") domain of red dots. 




Fig. 2.2. An example where the optimal integral solution (yellow dot) is far away from the 
optimal solution of the integer programming relaxation (white dot) and its nearest integer-feasible 
solution (red dot). 



Proof. If the optimal objective value z* of (IP) is achieved at the point x, then x 
is feasible for (IP), and hence it is also feasible for (LP). Therefore, z > c T x = z* . □ 

The consequence of the second effect is that it is much easier to solve the problem 
(LP) than (IP). 

We will shortly turn to studying how to solve LPs. Before we do so, let us 
remark that solving the LP relaxation and rounding the optimal values of the decision 
variables to the nearest feasible integer valued feasible solution is not always a good 
approximation scheme for integer programming problems, as the graphical example 
of Figure 2.2 shows. 

2.1. The Simplex Method Explained by Example. We will now discuss 
the simplex algorithm for solving general LPs. To motivate the method, we start 
with a simple example: 
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Example 9. Consider the LP instance 



z = max 5a?i + 4x2 + 3x3 

X 

s.t. 2x\ + 3x2 + xs < 5 
4xi + x 2 + 2x 3 < 11 
3xi + 4x 2 + 2x 3 < 8 
xi,x 2 ,x 3 > 0. 



Preliminary step I: introduce slack variables 

z = max 5xi + 4x2 + 3x 3 + OX4 + OX5 + 0x6 
s.t. 2xi + 3x2 + x 3 + X4 = 5 
4xi + x 2 + 2x 3 + x 5 = 11 
3xi + 4x 2 + 2x 3 + x 6 = 8 

Xl,X 2 ,X 3 ,X4,X 5 ,X 6 > 0. 

Preliminary step II: express in dictionary form 

max z s.t. xi , . . . , xq > 0, where 

X4 = 5 — 2xi — 3x 2 — x 3 
x 5 = 11 — 4xi — x 2 — 2x 3 
x 6 = 8 — 3xi — 4x 2 — 2x 3 
z = + 5xi + 4x 2 + 3x 3 . 

Step 0: xi,X2,x 3 = 0, x 4 = 5, x 5 = 11, x 6 = 8 is an initial feasible solution. 
xi,x 2 ,x 3 are called the nonbasic variables and X4,x 5 ,x 6 basic variables. That is, ba- 
sic variables are expressed in terms of nonbasic ones. 

Step 1: We note that as long as xi is increased by at most 

5 . /5 11 8 
- = mm - 



2 V 2 ' 4 ' 3 

all Xi remain nonnegative, but z increases. Setting x\ — 5/2 and substituting into the 
dictionary, we find X2,x 3 ,X4 = 0, X5 = 1, xq — 1/2, z — 25/2 as an improved feasible 
solution. We call xi the pivot of the iteration. 

We can now express the variables x\, X5, xq, z in terms of the new nonbasic vari- 
ables X2,x 3 ,X4 (those currently set to zero) to obtain a new dictionary. To do this, 
use line 1 (the line that relates xi and X4) of the dictionary to express xi in terms of 

X2,X 3 ,X 4 , 

xi = i (5 - 3x 2 - x 3 - x 4 ) 

and substitute the right hand side for x\ in the remaining equations. The new dictio- 
nary then looks as follows, 
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X\ 



5 3 1 

2 " 2 X2 2 
.x 5 = 1 + 5x 2 + 2x 4 



^3 



1 1 



1 



" 6= 2 + 2 2;2 "2" 34 
25 7 1 

Z= Y~2 X2 + 2 X3 



-X4 



2 Xi 



2 Xi - 



(2.1) 
(2.2) 
(2.3) 

(2.4) 



Of course, we are still solving 



max z s.t. x\, . . . ,xq > 0, 



subject to the relationships (2.1)-(2.4) holding between the variables, and the new 
LP instance is equivalent to the old one. However, a better feasible solution can be 
read off the new dictionary by setting the nonbasic variables to zero! 



Step 2: We continue in the same vein: increasing the value of xi or X4 is useless, as 
this would decrease the objective value z. Thus, £3 is our pivot, and we can increase 
its value up to 

1 = min(5, +00, 1), 

leading to the improved solution xi^x^x§ = 0, Xi = 2, £3 = 1, X5 = 1, z = 13 and 
the dictionary corresponding to X2,X4,xq as nonbasic variables, 

£3 = 1 + x 2 + 3x 4 — 2xq 
xi = 2 — 2x 2 — 2x 4 + x 6 

X5 = 1 + 5X2 + 2x4 
z = 13 — 3X2 — X4 — Xq. 



At this point we can stop the algorithm for the following reasons: From the last 
line of the dictionary we see that for any strictly positive value of X2,X4 or xq the 
objective value z is necessarily strictly smaller than 13, and from the other lines of the 
dictionary we see that as soon as the values of X2 , X4 and xg are fixed, the values of 
X3, x\ and x 5 are fixed too. Thus, the last dictionary yields a certificate of optimality 
for the identified solution. 



2.2. General Form of the Simplex Algorithm. Let us now try to understand 
how the dictionary 

£3 = 1 + x 2 + 3x4 - 2x 6 

xi = 2 - 2a; 2 - 2a;4 + x 6 (2.5) 
£5 = 1 + 5x 2 + 2x 4 
z = 13 — 3^2 — X4 — x 6 , 
15 



which was obtained after two pivoting steps, could have been obtained directly from 
the input data of the original LP instance 

(LPI) max 5x i + 4x 2 + 3.x 3 + 0x 4 + 0x 5 + Oxq 
s.t. 2x\ + 3x2 + x 3 + X4 = 5 
4xi +x 2 + 2x 3 + x 5 = 11 
3xi + 4x 2 + 2x 3 + x 6 = 8 

Xi,X2,X 3 ,X4,X 5 ,X 6 > 

if we had identified the relevant basic variables: 

The constraints of (LPI) imply a functional dependence between the nonnegative 
decision variables Xi , expressed by the linear system 

Ax = b, (2.6) 

where 





"2 


3 


1 


1 





0" 




"5 


A = 


4 


1 


2 





1 





, b = 


11 




3 


4 


2 








1 
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The basic variables of dictionary (2.5) are £3, £1,2:5. Writing 
x B := [x 3 xi x 5 ] T , x N := [x 2 x 4 x 6 ] T 

for the vectors of basic and nonbasic variables respectively, and extracting the sub- 
matrices 





"1 


2 


0" 




"3 


1 


0" 


A B = 


2 


4 


1 


, A N = 


1 










2 


3 







4 





1 



consisting of the columns of A that correspond to basic and nonbasic variables in the 
same order, the linear system (2.6) can be written as 

A B x B + A N x N = b. 

Solving for the basic variables, we obtain 

x B = Ag 1 {b - A N x N ) . (2.7) 

Likewise, the objective function can be written as 

Z = CgX B + cjj X N , 

where 

c B =[3 5 0] T , cat =[4 0] T , 

and using (2.7), this can be expressed as 

z = CgAgH + (cjf - c^A^An) x n . 
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The dictionary (2.5) is therefore the same as the system of equations 

x B = A B l b - A B l A N x N , 
z = CgA B 1 b + (c N - c B A B 1 A N ) x n . 

More generally, let us consider a general LP instance 

(P) max c T x 

s.t. Ax = b, 

x>0, 

where A is a m x n matrix with m < n. 
Definition 2.3. 

i) A dictionary of (P) is a system of equations 

x B = A B x b - Ag X A N x N , 
z = CgA^b + (cjf - CbA^An) x n , 

equivalent to 

Ax = b, 

T 

z — c x, 

where xb is a subvector (the so called basis) of size m of x such that themxm 
matrix A b consisting of the columns of A that correspond to the components 
of xb is nonsingular, and where xn is the remaining part of X and An the 
submatrix of A consisting of the columns corresponding to Xn- 
ii) Furthermore, this is called a feasible dictionary if A^b > 0, so that x = 
(xb,xn) — (Ag^b, 0) is a feasible (but generally suboptimal) solution for (P). 
(xb,xn) is then called a basic feasible solution. 

Having formally defined the notion of dictionary, we are ready to describe the 
dictionary form of the simplex algorithm in full generality: 

Algorithm 2.4. 

1. Choose a basic feasible solution (xb,Xn). 

2. Until Cn — AjjAg T CB <0 (componentwise), repeat: 

i) Choose i e {£ e N : q > AJA^cb}, where Ag is the i-th column of A. 

ii) If Ag Ai < (componentwise) 

(P) is unbounded (objective — > oo), stop. 
Else choose j G {£ e B : 5/J 'fie = mm{Sk/fik ■ fJ-k > 0}}, where 5k is the 
k-th entry of A^b and fik the k-th entry of A~^ A- L . 
End. 

Hi) Move i from N to B and j from B to N . 
End 

3. (xb,xjm) — (Af^b, 0) is an optimal basic solution, stop. 
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Before proceeding further, a few remarks are in order: 

The algorithm may get stuck in an infinite loop in which the same sets {B, N) 
reappear after periodic cycles. One can prove that this can be avoided by always 
choosing the smallest available candidate indices for i and j in steps i) and ii) (Bland's 
Rule). 

In an efficient implementation the full dictionary is never computed, as it suffices 
to compute the vectors cat— AjfAg T CB and A~j^Ai (which computationally amounts to 
solving two linear systems, one matrix- vector multiplication and one vector addition), 
as well as A~^b (a clever implementation does not require an extra linear systems 
solve). 

Any LP instance can be formulated in the form (P) by adding slack variables 
to linear inequalities and by replacing unconstrained variables Xi by the difference 
Xi = Si — ti of two bound constrained variables Sj, U > 0. 

In cases where an initial basic feasible solution is not available one can first solve 
the auxiliary LP instance 

m 

(AUX) max -^(si+ti) 

y=(s,t,x)m m + m + n ^ 



s.t. [I -I A] 

s,t,x > 0. 



A basic feasible solution for (AUX) is readily given by 

B = {i : 1 < i < m, b, L > 0} U {to + i : 1 < i < m, b t < 0}, 

so that ys = \b\ consists of components of s and t. Further, the optimal basic 
solution y* = (s*,t*,x*) of the auxiliary problem satisfies exactly one of the following 
two conditions: 

i) Either s*,t* — and then w.l.o.g. all the y-indices corresponding to s and t 
can be assumed to have been pivoted into the set of nonbasic variables y* Nl so that 

(x B ,x N ) = {^y B ^xV*N) 

is a basic feasible solution for (P), where the maps n x are the projections of the 
relevant parts of y onto their x-componcnts. 

ii) Or else, X^i( s i + U) > proves that the set of feasible solutions for (P) is 
empty, since any feasible solution x of (P) together with s,t = would improve on 
the optimal objective value. 



2.3. The Tableau Format of the Simplex Method. An alternative to using 
dictionaries to represent the data in intermediate steps of the simplex algorithm are 
tableaux. For example, the initial dictionary of the LP instance considered in the 
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previous section, 



x 4 = 5 — 2xi — 3x 2 — X3 
x$ = 11 — Ax 1 — X2 — 2x3 
^6 = 8 — 3xi — 4^2 — 2x3 
z = + 5xi + 4x2 + 3x3, 

can be written in tableau format as follows, 



2 


3 


1 


1 








5 


4 


1 


2 





1 





11 


3 


4 


2 








1 


8 


5 


4 


3 











0. 



The tableau is obtained by moving all the Xj in the dictionary to the left-hand 
side and constants to the right-hand side, then multiplying the last row by — 1 and 
omitting all the variables (and in the case of z, its coefficient —1). 

Note that the basic variables can be identified via an identity matrix that appears 
in the tableau. Furthermore, the rightmost entry of the last row equals the negative 
of the objective value at the current basic feasible dictionary. Obviously, an optimal 
dictionary corresponds to a tableau with all nonpositive entries on the last row (apart 
from possibly the rightmost entry). 

Applying the same procedure to the second dictionary 



5 3 1 1 

x 5 = 1 + 5x 2 + 2x 4 
1 

2 ' 2~" 2 

25 7 1 5 

y- 2*2 +2*3 -3*4, 



x e 



1 1 3 

X2 - 772:3 + 77*4 



we obtain the following tableau, 



1 


1.5 


0.5 


0.5 








2.5 





-5 





-2 


1 





1 





-0.5 


0.5 


-2.5 





1 


0.5 





-3.5 


0.5 


-2.5 








-12.5. 



Here is how the first tableau can be directly transformed into the second one (it is 
straightforward to check that this is the tableau equivalent to chaning the dictionaries 
according to the rules of the simplex algorithm) : 

1. Among the columns on the left, identify one whose last entry is positive. 
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2 


3 1 


1 


o 


o 


5 


4 


1 2 





1 





11 


3 


4 2 








1 


8 


5 


4 3 











0. 



We call this column the pivot column. 

2. For each row (apart from the last) for which the entry t in the pivot column 
is positive, look up the entry u in the rightmost column. If no such row exists, the 
problem is unbounded. If several such rows exist, pick the one for which t/u is small- 
est and call it the pivot row. 



2 


3 1 


1 








5 


4 


1 2 





1 





11 


3 


4 2 








1 


8 


5 


4 3 











0. 



3. Divide the pivot row by t, 



1 


1.5 


0.5 


0.5 








2.5 


4 


1 


2 





1 





11 


3 


4 


2 








1 


8 


5 


4 


3 











0. 



4. For all other rows i of the tableau, subtract the "rescaled" pivot row t; t times, 
where U is the row-z entry of the pivot colum, 



1 


1.5 


0.5 


0.5 








2.5 





-5 





-2 


1 





1 





-0.5 


0.5 


-2.5 





1 


0.5 





-3.5 


0.5 


-2.5 








-12.5. 



2.4. First Graphic Interpretation of the Simplex Method. For the pur- 
poses of understanding how the simplex method works, it is helpful to interpret it 
geometrically. We give two such interpretation, the first being more natural, while 
the second is the historical precedent and gave the method its name. 

Consider the polyhedron 

P = {x £ R n : Ax = b, x > 0} 
that constitutes the feasible domain of the LP problem 

(P) max c T x 

s.t. Ax = b, 
x>0. 
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Fig. 2.3. First graphical interpretation of the simplex method: Iterates move greedily from 
vertex to vertex so as to increase the objective function for as long as possible. 



Each set of nonbasic variables xn of a basic feasible solution identifies a vertex (ex- 
treme point) of P via the system of equations 



~A N A B 




x N 




b 


I 




X B _ 








Each pivoting operation changes N by only one entry and moves to a neighbouring 
vertex where the objective value is greater. The move can be seen as occurring along 
the ID facet 

{iel": Ax = b, x e = W G N \ {i}}. 

This graphic interpretation of the simplex method illustrates that basic feasible 
solutions correspond to vertices of the feasible polyhedron. 

It is trivial to see geometrically that generically an LP has a unique optimal so- 
lution located at a vertex, and that if the optimal solution is not unique (this can 
happen when a facet of the feasible polyhedron is orthogonal to the objective vector 
c), then the set of optimal solutions contains a vertex. 

Corollary 2.5 (Fundamental Theorem of Linear Programming). An LP has 
an optimal solution x* if and only if it has an optimal basic feasible solution. 

2.5. Second Graphic Interpretation of the Simplex Method. For the 

second graphic interpretation we consider an LP instance with an additional convexity 
constraint 

(CP) max c T x 

s.t. Ax = b, 
l T x = 1, 

x > 0, 
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Fig. 2.4. Second graphical interpretation of the simplex method: Simplices are "tumbling" up. 

where 1 = [ 1 ... 1 ] T . 

Writing z = c T x, we define (m + l)-dimensional points 





'b 




'A{ 




A n 


V = 






,...,w„ = 




z 







where Ai is the i-th column of A as before. The objective of (CP) is thus to express 
v as a convex combination of Wi, . . . , w n so as to maximise z. 

Due to the convexity constraint l T x = 1, a basic feasible solution is associated 
with a collection of m points {wj : j E B} for which {Aj : j £ B} are linearly 
independent and contain b in their convex hull. This is the same as requiring that 
the line {[^] : ( g 1} intersects the simplex conv({wj : j E B}). Thus, the simplex 
algorithm proceeds via a sequence of simplices that are "tumbling upwards" , and in 
each pivoting operation only one vertex of the simplex changes. 

2.6. LP Duality. Let us again consider the LP instance we studied previously, 

(P) max 5xi + 4^2 + 3x3 

s.t. 2x\ + 3x2 + %3 < 5 
4xi +x 2 + 2x 3 < 11 
3xi + 4x 2 + 2x 3 < 8 
xi,x 2 ,x 3 > 0. 

We saw that the optimal value is 13. 

In integer programming, instead of solving an LP relaxation to optimality one is 
often interested in finding merely upper and lower bounds on the optimal value. A 
lower bound is provided by any feasible solution. For example, xi, X2 = 1, X3 = is 
feasible with objective value 9. 

How can we obtain upper bounds? Multiplying the first constraint by 3 we obtain 

6x1 + 9x2 + 3x3 < 15, 
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and since Xi,X2,x 3 > 0, this yields an upper bound on the objective function: 
z = 5x\ + 4x 2 + 3x 3 < 6x1 + Qx 2 + 3x 3 < 15, 

Likewise, taking the sum of the first two constraints yields the valid upper bound 
z = hxi + 4:X 2 + 3x 3 < 6x1 + 4x 2 + < 16. 

More generally, such bounds can be obtained from any sum of positive multi- 
ples of the constraints for which the resulting coefficients are no smaller than the 
corresponding coefficients of the objective function: 



[5 4 3] < [yi V 2 2/3] 



2 3 1 
4 1 2 

3 4 2 



, Vi, 2/2, 2/3 > => z < [2/1 2/2 2/3] 



The best such upper bound is obtained by solving the LP instance 

(D) min52/i + H2/2 + 8y 3 

v 

s.t. 2 Vl + 4y 2 + 3y 3 > 5, 
32/1 + 2/2 + 42/3 > 4, 
2/i + 2y 2 + 2y 3 > 3, 
2/1,2/2,2/3 > 0. 

Definition 2.6. (D) is called the dual of (P), while (P) is called the primal 
problem. 

Using analogous arguments, it is easily worked that that a general LP of the form 

(P) max c^x + cjs 

(x,s) 

s.t. j4^x ~\~ -^4^^ = 
B x x + B s s < 6, 
x > 



is associated with a dual 

(D) min a T y + b T t 

(y,t) 

s.t. Ajy + B^t > c x 
A T s y + Bjt = c s 
t > 0. 

Furthermore, the dual of (D) is the primal (P). 

Let us now work out the relationship between primal and dual dictionaries. This 
will lead to an exact characterisation of optimal solutions of LPs. We consider the 
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primal-dual pair of LPs in general form, 

(P) max z = c T x O (D) min w = b T y 

x y 

s.t. Ax < b, s.t. A T y > c, 

x>0, y>0. 

Note (D) is equivalent to the maximisation problem 

(D') max — w = —b T y 
v 

s.t. - A T y < -c, 

y>0 

so that we can apply the simplex algorithm in the form we discussed. 

Let us call the variables of (P) x m +i, . . . , x m+n and introduce slack variables 
a;i,..., x m . For (D) on the other hand we introduce slack variables y m +i, . . . , y m +n 
and call the regular variables y\, . . . , y m . 

Theorem 2.7. Each dictionary 

x r = b r — UrsXs, (r e B) 

z = d+^2 Cs%s 

of the primal problem with augmented decision vector x, 

(P') max z — c T x 

X 

s.t. [I A]x = b, 

x>0 



corresponds to a dictionary 



—w = 



-C s + ^2 a rsVr, (s G N) 
r<=B 



reB 

of the corresponding augmented counterpart of (D 7 ), 

(D") max — w = —b T y 
v 

s.t. [-A T I]y = -c, 

y>o. 
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Let us illustrate Theorem 2.7 with an example: The initial dictionary 



Xi = 5 — 2^4 — 3x 5 — x 6 
X2 = 11 — 4x 4 — x 5 — 2x 6 
x 3 = 8 — 3x4 — 4x 5 — 2x 6 
z = + 5x 4 + 4x 5 + 3x 6 

of the problem 



(P') max 5x 4 + 4xs + 3x6 

s.t. 2x 4 + 3x5 + x 6 + xi = 
4x 4 + x 5 + 2x 6 + x 2 = 
3x 4 + 4x 5 + 2x 6 + x 3 
x 4 , . . . , x 6 > 

considered earlier corresponds to the dual dictionary 

y 4 = -5 + 2j/i + 4y 2 + 3y 3 
2/5 = -4 + 3yi +V2+ 4j/ 3 
2/6 = -3 + I2/1 + 2y 2 + 2y 3 
-to = - 5j/i - llj/2 - 8y 3 

of 

(D") max - 5yi - llj/2 - 8y 3 

y 

s.t. - 2j/i - 4y 2 - 3y 3 + 2/4 = -5, 

- 32/i - 2/2 - 4y 3 + j/5 = -4, 

- 2/i - 2y 2 - 2y 3 + 2/6 = -3, 
2/i,..., 2/6 > 0. 

Note that the dual dictionary is not feasible, as it corresponds to negative values for 
2/4, • • • , 2/6- 

The situation is very different with the final primal dictionary, 

x 6 = 1 + x 5 + 3xi — 2x 3 
x 4 = 2 — 2x 5 — 2xi + x 3 
x 2 = 1 + 5x5 + 2xi 
z = 13 — 3x 5 — xi — x 3 , 

which corresponds to the dual dictionary 

2/5 = 3- ?y 6 + 2j/ 4 - 5y 2 
2/1 = 1 - 3?/6 + 2y 4 - 2y 2 
2/3 = 1 + 2y 6 - 2/4 
-w = -13 - 2/6 - 2y 4 - 2/2- 



= 5 
= 11 

= 8 



This dictionary is not only dual feasible, but also optimal. It is immediate to see that 
this would always be the c < and b > 0. 

This leads to the following corollaries: 

Theorem 2.8 (LP Duality Theorem). If (P) has an optimal basic feasible so- 
lution x* then (D) has an optimal basic feasible solution y* too, and the primal and 
dual objective values coincide, 

T * iT * 

c x —by. 



Theorem 2.9. [Complementary Slackness] A feasible solution x\, . . . ,x* of (P) 
is optimal if and only if there exists a dual (optimal) feasible solution y*,..., y^ such 
that 

m 

x* > => ajjy* = c 3 

i=l 

n 

aijX* < bj y* = 0. 

Complementary slackness is a useful tool to check whether a given primal feasible 
solution is optimal. 

Example 10. Consider the LP instance 

(P) max5a;i + Ax 2 + 3x 3 

s.t. 2xi + 3x2 + x% < 5 
4xi +x 2 + 2x 3 < 11 
3xi + 4x 2 + 2x 3 < 8 

X!,X 2 ,X 3 > 0. 

Assuming that we somehow suspect that the solution 2 = 0, Xg = 1 is optimal 

for (P), here is how we can verify this claim: 

First we check that x* is feasible for (P) by substitution into the constraint inequal- 
ities. Second, if x* is optimal, it should be possible to apply Theorem 2.9to construct 
a dual feasible solution y* such that b T y* = c T x* . 

By complementary slackness, x* , x^ > implies 

2y{ + 4y* 2 + 3y 3 * = 5 (2.8) 
y{ + 2y* 2 + 2y* 3 = 3. (2.9) 

Furthermore, 4x* + x* 2 + 2xg = 10 < 11 implies y\ = 0. Substituting into (2.8) and 
(2.9), we obtain 

tyi + 3^ = 5 
2/i* + 2y 3 * = 3 
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Solving this linear system, we find y\ = 1, = 1. By construction, y* satisfies the 
constraints A T y* > c, and we also see that y* > 0. Furthermore, 

c T x* = 5- 2 + 4- + 3- l = 13 = 5- l + ll-0 + 8- l = b T y*. 

This confirms the optimality of both x* and y* . 

3. Branch and Bound. Having studied LP relaxations of IPs and how to solve 
them using the simplex algorithm, we are now ready to apply these techniques in a 
framework designed to solve completely general integer programming problems: LP 
based branch-and-bound. 

The main idea of branch and bound consists in dividing a computationally hard 
problem into easier subproblcms and systematically exploit the information gained 
from solving these subproblcms. It is thus a divide and conquer approach. 

PROPOSITION 3.1. Consider the problem 

z = max{c T :r : x G 

where & denotes the set of feasible solutions (as we saw, & is usually defined im- 
plicitly via constraints). If & can be decomposed into a union of simpler sets & = 
^lU-'U^i and if 

z \j] .— m ax{c T a; : x e <^j} (j = 1, . . . , k), 

then 

z = maxz^. 

j 

Let us work through an example to see how such decompositions might arise in 
practice: 

Example 11. Let & be the set of feasible tours of the travelling salesman problem 
on a network of 4 cities. Let node 1 be the departure city. 

& can be subdivided & = =^"(1,2) U ^(1,3) U ^(1.4) into the disjoint sets of tours 
that start with an arc (1,2), (1,3) or (1,4) respectively. 

Each of the sets ^(1,2), ^(1.3) an d =^"(1,4) can be further subdivided according to 
the choice of the second arc, ^"(1,2) = =^"(i,2)(2,3) U ^(i,2)(2,4) etc. 

Finally, we see that each of these sets corresponds to a specific TSP tour and 
cannot be further subdivided. We have found an enumeration tree of the TSP tours. 

We see that & was decomposed on a first level, and then each of the constituent 
parts was further decomposed on a further level and so on. Also, each corresponds 
to a branch (or subtree) of the enumeration tree. 

Thus, Proposition 3.1 allows us to decompose a hard problem into a possibly large 
number of easier branch problems, and to find an optimal solution of & by comparing 
the solutions found for the branch problems. 
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Fig. 3.1. The implicit enumeration tree of Example 11. 



Of course, for even quite moderately sized problems such a tree can no longer be 
explicitly enumerated, as the number of leaves grows exponentially in the problem 
size. The idea of implicit enumeration is based on building up the enumeration tree 
as we explore it, and to prune certain parts that are not worth looking at before those 
parts are even generated. The pruning mechanisms are based on the following insight: 

PROPOSITION 3.2. Consider the problem 

z = max{c T x : x € J?}, 

and let = &\ U • • • U J?fc be a decomposition of its feasible domain into smaller sets. 
Let 2^1 < z \i\ < ^b'l be lower and upper bounds on z^l = max{c T x : x G J?j} for all 
j. Then 

z := max z^ < z < max z^ =: z 
j j 

gives an upper and lower bound on z. 

Proof. For all j, z => z^ > z^\ so z > maxj zj^. Furthermore, z = max., z^ < 
maxjZ^l. □ 

3.1. Pruning Mechanisms. We speak of "pruning" a branch when we detect 
that we need no longer explore it further. This can happen for a variety of reasons. 

1. Pruning by Bound: A branch ^ can be pruned when z^ < z. 

2. Pruning by Infeasibility: If ^ = 0, then the corresponding branch can be 
pruned. 

Why would an empty subset &j appear in the decomposition of 
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Fig. 3.2. Pruning by bound. 




Fig. 3.3. Pruning by optimality. 

Remember that we don't set up the decomposition tree explicitly but use an 
implicit enumeration, typically by introducing more and more constraints as we move 
towards the leaves of the enumeration tree. 

As we proceed, the constraints may become incompatible and correspond to an 
empty set JPj, but this may not be a-priori obvious and has to be detected. For 
example, in the case where the bounds z^ , are computed via LP relaxation, the 
simplex algorithm will detect if the LP relaxation is infeasible, and then &j is empty 
as well. 

3. Pruning by Optimality: When z^ = z^ for some j, then the branch corre- 
sponding to needs no further consideration, since an optimal solution z^ = z^ = 
z^ for this branch is already available. Of course, we will not throw this solution 
away, as it may later turn out to be optimal for the parent problem &\ 

3.2. The Branch and Bound Concept. A branch and bound system system- 
atically exploits all of the above ideas to break down the solution of a hard optimisation 
problem max{c T x : x e into easier parts. Its main features are the following: 
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Implicitly build an enumeration tree for & by subdividing & into subproblems 
J^i U • • • U ^ and recursively proceeding the same with every &j . Build each part of 
the tree only when it is actually being explored. 

For each subproblem - corresponding to a branch of the enumeration tree 
- we compute dual (and possibly primal) bounds on the optimal objective value. 
Relaxation by tractable problems is used to compute dual bounds, while heuristics 
are used to compute primal bounds. 

Use Proposition 3.2 to tighten bounds at the root, and prune suboptimal branches 
before they are fully built. In other words, only build the parts of the enumeration 
tree that are actually worth exploring. 

If LP relaxation is used for generating dual bounds in a branch-and-bound sys- 
tem, we speak of LP based branch-and-bound. We will now explore this framework 
in more detail: Throughout the algorithm we will maintain and update a list AN of 
active nodes (denoted by their feasible sets we denote the feasible sets of their 
LP relaxations by Pj), a primal bound z on max{c T x : x G and an incumbent 
x* (the best solution encountered so far). 

3.3. An Example of Branch and Bound in Action. As a way of illustrating 
the LP based branch-and-bound technique, let us now solve the IP instance 



of . Introducing slack variables X3 , X4 , X5 and applying the simplex algorithm, we 
obtain the optimal dictionary 



z = max4xi — x 2 

s.t. 7xi - 2x 2 < 14 
x 2 < 3 

2xi — 2x2 < 3 
x G TL\. 



1. Bounding: To obtain a first upper bound, we solve the LP relaxation 



(P) z = max4xi — x 2 

s.t. 7xi - 2x 2 < 14 



x<i < 3 

2xi - 2x 2 < 3 
x > 



20 1 



2 





x 2 = 3 - • x 3 
23 2 




from which we read off the nonintcgral solution (xi, x 2 ) = (^,3) that corresponds to 
the upper (or dual) bound z = -y . 
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We could have generated a lower bound by finding a feasible solution for (IP) via 
a heuristic, but since we did not, we set z = — oo. 

2. Fractional branching: Since z < z, (J?) is not solved to optimality yet, so we 
need to branch. Since x\ is fractional, we distinguish the cases where x\ < |_xij and 
x\ > \x{\ . More generally, we can pick any index j such that Xj ^ Z and set 

J^i = J 5 " n {x : xj < [xj]}, 
= & n {x : xj > \xj~\}. 

Of course, futher down the tree we can use the same branching rule to subdivide any 
&j thus generated. In our case this leads to the two subproblems 

(J^i) z = max4xi — X2 

s.t. 7xi - 2x 2 < 14 
x 2 < 3 

2xi - 2x 2 < 3 
xi < 2 

x e z^. 

and 

(=^2) 2 = max4xi — x 2 

s.t. 7xi - 2x 2 < 14 
x 2 < 3 

2xi - 2x 2 < 3 

xi > 3 

x e z^. 

Note that upon solving the LP relaxations of these two problems we will have 

max{z [1] , z [2] } < z, 

since x\ would have to be allowed to take the value ^ for the value ~z to be attained. 
Thus, branching on a fractional variable guarantees that the upper bounds decrease 
strictly in subsequent iterations, unless there exist multiple optimal solutions for the 
parent LP relaxation. 

So far we obtain the partial enumeration tree of Figure 3.4. The branches dangling 
off the nodes &\ and ^2 still need to be explored. We mark these nodes as active. 
The node & on the other hand has been processes and is inactive. 

3. Choosing an an active node for processing: During the run of the algorithm, 
we maintain a list of active nodes. This list currently consists of ^1,^2- Later we will 
discuss breadth-first versus depth-first choices of the next active node to be processed. 
For now we arbitrarily select &\ . 
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xi < 2 



xi > 3 



Fig. 3.4. The partial enumeration tree after the first fractional branching. 

4- Bounding: Next we derive a bound z' 1 ' by solving the LP relaxation 

(Pi) z = max 4a; i — X2 

s.t. 7xx - 2x 2 < 14 
x 2 < 3 

2 Xl - 2x 2 < 3 
xi < 2 
x>0 

of problem = max{c T a; : x G 

Note that (Pi) differs from (P) by having one more constraint x\ < 2 imposed. 
In the next section we will see that such problems can be solved via a warm start 
procedure that works out as follows: Using x\ = ^ — jx 3 — jX4 from the optimal 
primal dictionary of (P), express the constraint x\ < 2 as 

20 1 2 

— - lfX 3 - -Xi < 2. 

Introduce a new slack variable xq and add the constraint to the dictionary. 
The old dual solution is still dual feasible for the new dictionary 





20 


1 2 


Xl = 


y 


- ^xz - ^x 4 


x 2 = 


3 - 


• X3 — X4 




23 


2 10 


x 5 = 


y 


+ -jX Z - —x A 




6 


1 2 


x 6 = 


~7 


+ + 




59 


4 1 


z = 


y 


- y.X 3 - y£4- 



After just a few pivots on the dual dictionary and translating the optimal dual 
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xi < 2 



x\ > 3 




Fig. 3.5. The partial enumeration tree after the second fractional branching. 



dictionary back into a primal dictionary, we obtain 



Xi 


= 2-x 6 




1 1 


x-i 


= j + 2 X5 + X6 


x 3 


= 1 + x 5 + 5x 6 




5 j. 1 ^fi 


x 4 


= 2 + 2^5 + 6a; 6 




= y - ^ _ 3x6 



Therefore, zM = f , (af 1 ,^ 1 ) = (2, §). 

5. Branching: (x^ , a;^) is still not integral. Hence, J^i is not solved to optimality 
and cannot be pruned. Instead, we need to branch further. 

Using the same approach as before, we introduce 

=^11 =^in{i: i2<0} = ^n{i: x 2 = 0} (since x 2 > 0), 
J? 12 =^fl{i: x 2 > 1}. 

We now arrive at the partial enumeration tree of Figure 3.5, and the new list of active 
nodes is ^11,^12,^2- 

6. Choosing an active node: We arbitrarily choose J^2 for processing. 
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7. Bounding: We compute a bound z' 2 ' by solving the LP relaxation 

(P 2 ) z = max4xi — x 2 

s.t. 7xi - 2x 2 < 14 
x 2 < 3 

2xi - 2x 2 < 3 
xi > 3 
x > 



of the problem 



(^2) z [1] = max{c T x : x G J^}- 



Proceeding the same way as above, we express the new constraint x\ > 3 in terms 
of the optimal nonbasic variables of (LP). We had found xi = 20/7 — 1/7x3 — 2/7x 4 . 
Thus, the new constraint is 

1 2 1 

-x 3 + -x 4 < --. 

After introducing a new slack variable Xq, we find that (LP 2 ) is equivalent to 

59 4 1 
z = max — — -x 3 — -Xi 

20 1 2 

S.t. Xl = y - -X 3 - -X 4 

X2 = 3 — • X3 — Xi 
23 2 10 

£5 = — + jX 3 - —Xi 

1 1 2 
xi, . . . ,x 6 > 0. 

But this LP is infeasible, because the last constraint contradicts x 6 > 0. Hence, & 2 
can be pruned by infeasibility, and we arrive at the partial enumeration tree of Figure 
3.6. 

8. Choosing an active node: From the list {^11,^12} of active nodes we arbi- 
trarily choose ^12. 

9. Bounding: Proceeding as above, we solve the LP relaxation of (^12) and find 
the optimal solution (xf 2 ' , x 2 12 ' ) = (2, 1) corresponding to the optimal value z' 12 ' = 7. 
Since x' 12 ' is integral, this is a feasible solution for (J^i) and provides a lower bound 
zj 12 l = 7. In fact, ^i 2 can now be pruned by optimality. The partial enumeration 
tree is now as in Figure 3.7. 

10. Updating the incumbent: We store (1, 2) as the best integer solution found so 
far and update the lower bounds z <— max(z, 7) = 7, zj 1 ^ <— maxfz' 1 ' , 7) = 7, 

11. Choosing an active node: Only is active, so choose this node. 

12. Bounding: Proceeding as above, we solve the LP relaxation of (^11) and find 
the optimal solution (x^ 11 ',^ 11 ') = (f ,0) with optimal value z' 11 ' = 6. 

34 




Since z' 11 ' < z, we can prune &\\ by bound and arrive at the partial enumeration 
tree of Figure 3.8. 

13. Termination: There are no further active nodes left, and the algorithm 
terminates, returning the optimal solution z = 7 and the maximiser x — (2, 1) that 
achieves it. 

3.4. LP Based Branch and Bound: The General Case. Let us now de- 
scribe the LP based branch-and-bound algorithm for general IP instances. We start 
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Fig. 3.8. Pruning j^ii by bound yields the terminal partial enumeration tree. 



with a variant that only uses dual bounds. 

Algorithm 3.3. [LP based branch- and-bound] 

AN = {.^}, z = — oo, x* — (initialisation) 
while AN 7^ 0, repeat 

choose a problem &j e AN and apply simplex alg. to x^ := argmax{c T a; 

x e Pj} 

if simplex alg. found Pj = 
set AN = AN\{&j} 

elseif simplex alg. found LP relaxation unbounded 

set = z = +oo and AN = (AN \ {&j}) U {&f\&f ] } (branch- 
ing) 

else (x^ is finite) 

set = c T x^\ z = max^I^) 
ifzW < z 

set AN = AN\{&j} 
elseif is feasible for & 

set z = x* = xW and AN = AN \ : < z} 
else (&j not pruned) 

set AN = (AN \ {J^}) U {&f\&f ] } (branching) 

end 

end 

end. 
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Remarks: 



The algorithm is constructed to either stop with an optimal incumbent x* , or 
with x* = as certificate of infeasibility of the problem. In the first case the final 
bounds satisfy z = ~z. 

The efficiency of the algorithm can be improved if we also compute a lower bound 
for &j (using a heuristic) and use this to update the global lower bound z. 

In LP-based branch-and-bound the LP relaxation of each subproblcm is only 
slightly different from that of its parent problem. In this situation the simplex algo- 
rithm does not have to be started from scratch. Instead, it is much more economical 
to use a warm start from a basic feasible solution constructed from the optimal dic- 
tionary of its parent. We will discuss this technique in the next section. 

Algorithm 3.4. [Variant with lower bounds] 

AN = {^}, z — ~oo, x* — (initialisation) 
while AN ^ 0, repeat 

choose a problem j£j e AN and apply simplex alg. to x 1 ^ := argmax{c T x : 

xePj} 

if simplex alg. found Pj = 
set AN = AN\{&j} 

else 

using heuristic, compute lower bound z^ 

set z = max{z,zW} and AN = AN \ : < z} 

if simplex alg. found LP relaxation unbounded 

set zW = z = +oo and AN = {AN \ {J? }) U {&f , &f ] } 
(branching) 
else (x 1 ^ is finite) 

set z^ = c T x^, z — max(l, z^ ) 
ifz® <z 

set AN = AN\ {J? } 
elseif x^ is feasible for & 

set z = zW , x * = xW and AN = AN \ {&, L : z^ < z} 
else (&j not pruned) 

set AN = (AN \ {^}) U {&f , &f ] } (branching) 
end 

end 

end 

end. 

4. Improved Branch and Bound Systems. Let us now discuss an array of 
tools that can be used to design more efficient branch-and-bound algorithms. 

4.1. Warm-Starting the LP Relaxations. In LP-based branch-and-bound 
the LP relaxation of each subproblem is only slightly different from that of its par- 
ent problem. In the example of Section 3.3 we have seen that in this situation the 
simplex algorithm does not have to be started from scratch. Instead, it is much more 
economical to use a warm start from a basic feasible solution constructed from the 
optimal dictionary of its parent. Let us now discuss this technique in more detail, 
distinguishing four different cases. Let x*,y* be optimal solutions for the primal and 
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dual LP instances 



(P) maxc T x (D) 

X 

s.t. Ax < b, 

x > 

4.1.1. The objective function changes. 

(P ! ) maxc T x 

X 

s.t. Ax < b, 

x > 0, 

Note that x* is still primal feasible. A few simplex iterations can therefore be 
applied to the primal dictionary corresponding to x* to find the new primal optimal 
solution x* . The new dual optimal solution y* can then be read off the new optimal 
dictionary. 

EXAMPLE 12. Let us once again consider the LP instance 

(P) max5xi + 4x 2 + 3x 3 

s.t. 2xi + 3x2 + x 3 < 5 
4xi + x 2 + 2x 3 < 11 
3xi + 4:X 2 + 2x 3 <8 
X1.x2.x3 > 

for which we already found the optimal dictionary 

x 3 = 1 + x 2 + 3x4 — 2x e 
x\ — 2 — 2x 2 — 2^4 + x e 
x 5 = 1 + 5^2 + 2x 4 
z = 13 — 3x 2 — x 4 — x 6 

corresponding to the primal optimal solution 3 = 1, x*2 = 0. Let us now 

change the objective function to z = 6x1 + 5x 2 + 2x 3 . This can be expressed in terms 
of the nonbasic variables 

z = 6(2 - 2x 2 - 2x 4 + x 6 ) + 5x 2 + 2(1 + x 2 + 3x 4 - 2x 6 ) 
= 14 — 5x 2 — 6x4 + 2x 6 . 

Thus, for the changed problem the dictionary corresponding to x\ = 2, X3 = 1, 
x* 2 = is as follows, 

x 3 = 1 + x 2 + 3x 4 - 2x6 
xi = 2 — 2x 2 — 2x 4 + x 6 
x 5 = 1 + 5x 2 + 2x 4 
z = 14 — 5x 2 — 6x4 + 2x6- 
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miri by 

y 

s.t. A T y > c, 

y>0- 



This dictionary is not optimal, as the nonbasic variable Xq appears with the positive 
coefficient 2. Thus we use x 6 as a pivot and increase it from zero to 

2 = mm(-,2,+oo). 

The new dictionary is 

11 13 

x& = g + 2 X2 ~ 2 X3 + 2 Xi 

5 3 1 1 
Xl = 2 " 2 s2 " 2 s3 " 2 X4 
x 5 = 1 + 5x 2 + 2x 4 
z = 15 — 4x 2 — x 3 — 3x 4 . 

77ms dictionary is optimal, as all the coefficients that appear in the z-line are nonpos- 
itive. We can read off the new primal optimal solution 



x l — 2' ^2 — x 3 — U ' 



and the new optimal dual solution 

Hi = 'i,V2 = yl = °- 

Both yield the new optimal objective function value 

6x{ + 5$2 + 25 3 = 15 = 5y{ + + 83/3. 



4.1.2. A new variable appears. 

(P') max c T x + c n+1 x n+1 

s.t. Ax + ax n+ i < b, 
x,x n+1 > 0, 

Note that (x*,0) is primal feasible for (P'). A few simplex iterations can therefore 
be applied to the primal dictionary corresponding to (x* , 0), and the new primal and 
dual optimal solutions can be read off the new optimal dictionary. 

Example 13. Consider the LP instance from Example 12 and let us add a new 
variable, 

(P') max5xi + 4x2 + 3x3 + X4 

s.t. 2x\ + 3x 2 + xs — X4 < 5 
4a; i + x 2 + 2x 3 + x 4 < 11 
3xi + 4x 2 + 2x 3 - 2x 4 < 8 
xi,x 2 ,x 3 ,x 4 > 0. 
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Introducing slack variables, 

£5 = 5 — 2x\ — 3x 2 — X3 + X4, 
^6 = 11 — 4a; 1 — x 2 — 2x3 — X4, 
x-j = 8 — 3.ti — 4x 2 — 2x 3 + 2x 4 . 

To find the dictionary corresponding to (a;*,0) = (2,0,1,0) (or (2,0,1,0,0,1,0) with 
the slack variables) we need to reformulate these equations so as to express x\,X3,Xe 
in terms of x 2 , X4, x 5 , X7. Gaussian elimination yields 

x\ = 2 - 2x 2 - 2x 5 + x-j, 
x 6 — 1 + 5x 2 — 3x4 + 2x 5 , 
x 3 = 1 + x 2 + X4 + 3x 5 - 2x 7 . 

Expressing z = 5x\ + 4x 2 + 3x 3 + x 4 in terms of x 2 ,x 4 ,X5,X7 by substituting 
from the previous equations, we find the required dictionary to warm-start the simplex 
method: 

x\ = 2 — 2x2 — 2x 5 + X7 
xg = 1 + 5x 2 — 3x4 + 2x5 
X3 = 1 + x 2 + X4 + 3x5 — 2x7 
z = 13 — 3x 2 + 4x4 — %5 ~ x-j. 

4.1.3. The right-hand side changes. 

(P') maxc T x 

X 

s.t. Ax < b, 

x > 0, 

Note that y* is still feasible for the new dual 

(D') min& T y 

y 

s.t. A T y > c, 

y>o, 

We can thus proceed as in Section 4.1.1 with (D') in the role of (P'). 

4.1.4. A new constraint appears. This is the most common change in the 
context of branching algorithms for IP. 

(P') maxc T x 

X 

s.t. Ax < b, 

a m+l x ^ &m+l) 

x > 0, 
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In the context of the new dual this corresponds to the situation where a new variable 
appears, 

(D') min b T y + b m+1 y m+ i 

y,y m +i 

s.t. A T y + a m+ iy m+1 > c, 
y,y m +i > 0. 

Thus, (y*,0) is feasible for (D'), and we can proceed as in Section 4.1.2 with (D') 
taking the place of (P')- 

4.2. Fine- Tuning Branch-and-Bound. Commercial branch-and-bound sys- 
tems for integer and mixed integer programming are designed to solve as few LP 
subproblems as possible and to solve them rapidly. Typically, such software allows 
the user to switch on or off a number of options to fine-tune the method to specific 
problem classes. We will now discuss a list of typical options and the mathematical 
reasoning behind them. 

4.2.1. Preprocessing. LP or IP models can often be simplified by reducing the 
number of variables and constraints, and IP models can be tightened before any ac- 
tual branch-and-bound computations are performed. 

Example 14. Consider the LP instance 

max 2xi + X2 — X3 
s.t. 5xi — 2x 2 + 8x 3 < 15 
8x1 + 3x 2 — x 3 > 9 
x\ + x 2 + £3 < 6 
< x x < 3 
< x 2 < 1 
l<x 3 . 

Tightening bounds: Isolating x\ in the first constraint and using X2 < 1, —^3 < — 1 
yields 

5xi < 15 + 2x 2 - 8x3 < 15 + 2 x 1 - 8 x 1 = 9, 

and hence, x\ < 9/5, which tightens the bound x\ < 3. 

Likewise, isolating X3 in the first constraint, and using the bound constraints, we 

find 

8x3 < 15 + 2x 2 - 5xi < 15 + 2 x 1 - 5 x = 17. 

This implies X3 < 17/8 and tightens X3 < 00. 

And finally, isolating X2 in the first constraint, 

2x 2 > 5xi + 8x3 - 15 > 5 x + 8 x 1 - 15 = -7 
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yields X2 > —7/2 which does not tighten X2 > 0. 

Proceeding similarly with the second and third constraints, we obtain the tight- 
ened bound 

8x1 > 9 - 3iE2 + x 3 > 9 - 3 + 1 = 7, 

yielding the improved bound x\ > 7/8. 

As some of the bounds have changed after the first sweep, we may now go back 
to the first constraint and tighten the bounds yet further. Isolating x 3 , we obtain 

7 101 

8x 3 < 15 + 2x 2 - 5a;i < 15 + 2 - 5 x - = — , 

8 8 

yielding the improved bound x 3 < 101/64. 

Continuing the second sweep by isolating each variable in turn in each of the con- 
straints 1-3, and using the bound constraints, several bound constraints may further 
tighten in general, but not in the present example. 

How many sweeps of this process are needed? One can show that after two sweeps 
of all the constraints and variables, the bounds cannot improve any further! 

Redundant Constraints: Using the final upper bounds in constraint 3, 

9 101 „ 

x 1 + x 2 + x 3 < - + 1 + — < 6, 
5 64 

so that this constraint is redundant and can be omitted. 
The remaining problem is 



max 2x i + x 2 — x 3 

5xi - 2x 2 + 8x 3 < 15 
8a; i + 3x 2 — x 3 > 



, _ 9 

7 9 101 

^<X!<-, 0<x 2 <l, l<x 3 <—. 

8 5 64 



Variable fixing: 

• Increasing X2 makes the objective function grow and loosens all constraints 
except X2 < 1. Therefore, in an optimal solution we must have X2 = 1. 

• Decreasing x 3 makes the objective function grow and loosens all constraints 
except 1 < x 3 . Thus, in an optimal solution we must have x 3 = 1. 

This leaves the trivial problem 



f 7 

max < 2xi : - < X\ < 



Example 14 shows how to simplify linear programming instances. In the prepro- 
cessing of IPs we have further possibilities: 

• For all Xj with an integrality constraint Xj £ Z any bounds lj < Xj < Uj can 
be tightened to < Xj < [uj\. 
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• For binary variables new logical or Boolean constraints can be derived that 
tighten the formulation and hence lead to fewer branching nodes in a branch- 
and-bound procedure. 

The latter point is illustrated in the next example: 

Example 15. Consider a BIP instance whose feasible set is defined by the fol- 
lowing constraints, 

7x\ + 3x2 — 4^3 — 2x4 < 1 
-2xi + 7x 2 + 3x 3 + x 4 < 6 
— 2x 2 — 3.t 3 — 6x4 < —5 
3xx - 2x 3 > -1 

x e B 4 . 

Generating logical inequalities: The first constraint shows that x\ — 1 =>■ £3 = 1, 
which can be written as 

x\ < x 3 . 

Likewise, x\ = 1 => X4 = 1, or equivalently, 

Xl < X4. 

Finally, constraint 1 also shows that the problem is infeasible if x\ = X2 = 1. 
Therefore, the following constraint must hold, 

Xl +x 2 < 1. 

We can process the remaining constraints in a similar vein: Constraint 2 yields the 
inequalities X2 < x\ and x 2 + x 3 < 1, constraint 3 yields X2 + X4 > 1 and x 3 + X4 > 1, 
and constraint 4 yields x\ > x 3 . 

Although the introduction of the new logical constraints makes the problem seem 
more complicated, the formulation becomes tighter and thus better. Furthermore, we 
can now process the problem further: 

Combining pairs of logical inequalities: We now consider pairs involving the same 
variables: The inequalities x\ < x 3 and x\ > x 3 yield xi = x 3l xi + x 2 < 1 and 
x 2 < x\ yield X2 = 0, and then X2 + X4 > 1 yields X4 = 1. 

Simplifying: Substituting the identities X2 — 0, x 3 — x\ and X4 = 1 we found, all 
four constraints become redundant. 

We are left with the choice x\ G {0, 1}, and hence the feasible set contains only 
two points 

S ={(1,0, 1,1), (0,0, 0,1)}. 

4.2.2. Selection of Algorithm to Solve the LPs. The simplex algorithm 
implemented in a branch-and-bound system typically lets the user choose some pa- 
rameters. Likewise, interior-point LP solvers can be fine-tuned via certain parameters. 
We cannot discuss this issue further in the framework of this course. 
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4.2.3. Node Selection. Depth-First: A depth-first strategy aims at finding a 
feasible solution quickly, by only choosing an active node that is a direct descendant 
of the previously processed node. This strategy also makes it easier to warm-start the 
LP calculations in each iteration. 

Best-Node-First: To minimise the total number of nodes processed during the run 
of the algorithm, the optimal strategy is to always choose the node with the largest 
upper bound, i.e., Sj such that 

z [ ' j] = max{z H : S { G List}. 

Under this strategy we will never branch on a node St whose upper bound z^ is 
smaller than the optimal value z of S. This is called a best-node-first strategy. 

The depth-first and best-node-first strategies are usually mutually contradictory, 
so a compromise has to be reached. Usually, depth-first is used initially until a feasible 
solution is found and a lower bound z is established. 

4.2.4. Branching Options. In our example, we branched on a fractional vari- 
able. Often several candidate variables are available. In this choice has to be 
made. 

Most Fractional Variable: Let C be the set of fractional variables of the solution 
x* to a LP relaxation. The most fractional variable approach is to branch on the 
variable that corresponds to the index 

j = argmaxmin{/j, 1 - /,}, 

iGC 

where ft = x* - [x*\. 

Branching by Priorities: In this approach the user can indicate a priority of im- 
portance for the decision variables to be integer. The system will then branch on the 
fractional variable with highest priority. 

Example 16. Consider a fixed charge network problem 

min{c T a; + fy: Nx = b, x < uy, x G M™ , y G Z™ }, 

where N is the node-arc incidence matrix of the network, and b is the demand vector. 

Since rounding up the variables yj corresponding to large fixed costs fj changes 
the objective function more severely than rounding yj corresponding to small fixed 
costs, we prioritise the yj in order of decreasing fixed costs fj. 

GUB/SOS Branching: Many IP models contain generalised upper bound (GUB) 
or special ordered sets (SOS) constraints of the form 

k 

with Xj G B for all j. In this case it is not good to branch on a factional variable x* 
from the solution x* of a LP relaxation, since the branching does not lead to balanced 
sets: 



Si = {x G S : Xi = 1} 
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contains only one point Xi = 1, Xj = Vj ^ i, whereas 

S 2 = {x e S : Xi = 0} 

contains |S| — 1 points. 

In these good choice of branching is given by fixing an order ji , . . . , jh of 

the variables and choosing 

51 =Sn{x: Xji =0, i = l,...,r}, 

5 2 = S n {x : x,-, = 0, i = r + 1, . . . , k}, 

where r = min{s : J2i=i ^ > 3}- 

Strong Branching: This is used only on difficult problems where it is worthwhile 
spending more time to find a good branching. In this approach one chooses a set C of 
candidate variables, branches up and down for each Xj <E C, computes upper bounds 
zf and z^ 3 by solving the LP relaxations of the up and down branching corresponding 
to xj , and chooses the variable having the largest effect 

j* = argmin{max(zj / , zf) : j E C} 

as the actual branching variable for the branch-and-bound scheme. That is, one 
explores several possible branchings and chooses to pursue the branches below only 
the most promising branching variable. 

4.2.5. Further Improvements to Branch-and-Bound Systems. LP based 
branch-and-bound may not be strong enough to solve really hard IP problems. In 
these cases one has to resort to the following improvements. 

Finding feasible solutions via heuristics. This is hard, but when it succeeds, it 
generates lower bounds that allow pruning by bound. 

Finding better dual bounds through combinatorial relaxation, duality, or semidef- 
inite programming relaxation. 

Tightening the formulation of the IP, having the effect that fewer branching nodes 
are needed until an optimal solution is found. This can be achieved with cutting plane 
and branch- and- cut methods, or by Lagrangian relaxation or column generation ap- 
proaches in which part of the feasible region is convexified. 

5. Alternative Formulations of Integer Programmes. Recall our model of 
alternative disjunctions from Section 1, 

T 

mm c x 

(x^)eR"+ 2 

s.t. ajx - b t < M(l - Vi ) (i = 1,2), 

Vi +2/2 = 1, 
< x < u 

y E B 2 . 
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It would have been mathematically equivalent to replace the constraint 

2/i + 2/2 = 1, 

corresponding to the requirement that exactly one of the alternative constraints ajx < 
bi (i = 1, 2) be imposed, by the constraint 

2/i + 2/2 > 1, 

corresponding to the requirement that at least one of the alternative constraints be 
imposed. 

Another example of an alternative formulation is the following: 

Example 17. Recall our binary programming model of the travelling salesman 
problem, 

n n 

min EE 

Cij Xij 

i=l 3=1 

s.t. ^ x ij = 1 {i € [l,n]), 
5^x tf =l (j'e[l,n]), 

x e B" x ™. 

TTie cut-set constraints 

were introduced to eliminate solutions that contain subtours. 

Alternatively, this could be achieved using subtour elimination constraints; 

EE^' - \ S \ ~ 1 V S C N s.t. 2 < \S\ < n — 1. 

ies jes 

The new model 

n n 

m i n E E ci o x a 

i=l j=l 

s.t. ^2 x v = 1 (* e i 1 ^])' 

^^iy<|S|-l V S C N s.t. 2 < \S\ < n — 1, 

ies jes 

x e M nxn . 
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Fig. 5.1. Subtours that don't visit all nodes need to be disallowed. 



is mathematically equivalent to the old model - the set of feasible tours has not 
changed, and neither has the objective function - but from an algorithmic point of 
view the formulation is very different. 

These examples raise the natural question as to whether it matters which for- 
mulation of an IP we use when several formulations are available. All integer pro- 
gramming problems have infinitely many different formulations that are all equivalent 
from a mathematical point of view in the sense that the optimal values of the decision 
variables are the same. However, the performance guarantee of algorithms can be 
very different for different formulations, so which formulation we use matters a great 
deal. 

Deriving alternative formulations of an IP is usually a nontrivial task. Conceptu- 
ally, many general approaches for solving IPs - including branch-and-bound - attempt 
to generate better and better formulations until a trivial formulation is found - one 
for which the optimal solution is obvious. 

5.1. Polytopes and Polyhedra. Before we discuss the business of alternative 
formulations further, we need to understand the geometric objects that are involved 
in defining them. 

Definition 5.1. 

i) A polyhedron is any subset o/K" described by a finite set of linear inequalities 

P = {x e R" : Ax < b}. 
ii) A polytope is the convex hull of finitely many points x\, . . . , Xk £ K n 



Hi) Let C C 1" be a convex set. A point x G C is an extreme point of C if x is 



k 




47 



Fig. 5.2. A polyhedron is the intersection of finitely many linear halfspaces. 




Fig. 5.3. A polytope is the convex hull of finitely many points. 

not a convex combination of two points in C distinct from x 

flxi,x 2 6 C\ {x}, A e (0, 1) s.t. x = Xxi + (1 - X)x 2 - 
iv) Let B C IR™. The convex hull of B is the set of points 

corw(B) = j^A^ : k e N, x t £ B, A, > 0, ^ A 4 = 1 j 

£/ia£ can 6e written as a convex combination of elements of B. 

The extreme points of a bounded convex set can be thought of as a minimal set 
of generators of the set because of the following theorem: 

Theorem 5.2 (Caratheodory's Theorem). Let C be a bounded convex subset of 
R n , and let E C C be the set of extreme points. Then the following are true: 

i) C = conv(E). 

ii) For every set D C C such that conv(D) — C we have E C D. 

Hi) Any point in C can be written as a convex combination of at most n + 1 
elements of E. 

It turns out that bounded polyhedrons are the same objects as polytopes. The 
characterisations in terms of constraints or generators are dual descriptions of one 
another 2 : 

Theorem 5.3. Let P be a bounded subset ofW 1 . Then P is a polytope if and 
only if it is a polyhedron. 

2 We note however that deriving one description explicitly from the other is a nontrivial compu- 
tational task 

48 



Fig. 5.4. A convex set and its extreme points (in red). 




Fig. 5.5. Polytopes can either be characterised by linear inequalities or by generators. 

5.2. Geometric Interpretation of Alternative Formulations. Let us now 

consider the mixed integer programming problem 

(MIP) max cEx + cTy 

s.t. Ax + By <b 

x G Z™, 

where A and B are matrices, and let us denote the set of feasible solutions of (MIP) 

by 

& := {(x, y) G M. n+P : Ax + By < b, x G Z"} . 

We call Z" x M p the integer lattice for this problem. 

Note that if we drop the integrality constraints x 6 Z n , the set of points that 
satisfy the remaining constraints is a polyhedron: 

P ■= {{x,y) e K" +p : Ax + By < b} , 

see Figure 5.6. The feasible set of the MIP is the intersection of this polyhedron with 
the integer lattice, & = P n (Z™ x M p ), see Figure 5.7. 

We now note that a reformulation of (MIP) generally defines a different polyhe- 
dron 

P := {(x, y) G M" +p : Ac + By < b} . 
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Fig. 5.6. The feasible set of the LP relaxation of an IP or MIP is a polyhedron. 



Fig. 5.7. The feasible set of the MIP is the intersection of the polyhedron with the integer lattice. 

The polyhedron P also has the property that its intersection with the integer lattice 
yields the feasible set of (MIP), & = Fn(Z" x W). In fact, there are infinitely many 
polyhedra with this property. But while the MIP described under an alternative for- 
mulation is the same, the LP relaxation is different, and algorithms such as LP based 
branch and bound behave differently. 

Definition 5.4. A formulation of the feasible set J? of the mixed integer pro- 
gramming problem (MIP) is any polyhedron P C R' i+P such that & — P(l (Z n x R p ) . 

Example 18. Consider the following set in B 4 , 

= {(0, 0, 0, 0), (1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1), (0, 1, 0, 1), (0, 0, 1, 1)} 

The following are different formulations of ^ , where 1 := [I, . . . , 1] T , 

Pi = {x e M 4 : < x < 1, 83xi + 6fx 2 + 49x 3 + 20x 4 < 100}, 
i 5 2 = {ieM 4 :0<i<l, 4xi + 3x 2 + 2x 3 + x 4 < 4}, 

F 3 = {i£l 4 : 4xi + 3x 2 + 2x 3 + x 4 < 4, xi + x 2 + x 3 < 1, xi + x 4 < 1, < x < 1} 

It can be seen that P 3 C Pi C P\ , that is, the formulations are successively tighter. 

Definition 5.5. Let & C M" be the feasible set of an IP or MIP, and let 
P,RC K™ be formulations of & ' . 

i) If P C R, we say that P is a tighter formulation than R. 
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Fig. 5.8. Two alternative formulations of the set of red dots. 




ii) If P — conv(^), we call P an ideal formulation. 

Tighter formulations of & are better, because algorithms tend to find an optimal 
solution faster. In particular, branch-and-bound systems are able to get away with 
far fewer branching steps. In the case of Example 18, it even turns out that P3 is an 
ideal formulation of & . In this case, the IP can be solved by solving the LP relaxation 
associated with the formulation P 3 . 



Theorem 5.6. If P C (K™ xl p ) is an ideal formulation of the mixed integer 
programming problem 

(MIP) z* = max cEx + cEy 

s.t. Ax + By <b 

x G Z'\ 

then any optimal basic feasible solution of the LP relaxation 

(LP) z LP = max cEx + cEy 
s.t. (x,y)eP 

is an optimal solution of (MIP). 

Proof. First, recall that any basic feasible solution (x, y) of (LP) corresponds to 
a vertex of P and thus is among the extremal points E of P. Since P = conv(j^"), 
Theorem 5.2 implies that x £ & , so that x is feasible for (MIP). Furthermore, (LP) 
being a relaxation of (MIP), we have zj,p = cjx + c^y > z* , so that (x,y) must in 
fact be optimal for (MIP). □ 
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6. Techniques for Finding Bounds on Integer Programmes. In our dis- 
cussion of the branch-and-bound approach we have learned that given an integer 
programming problem 

(IP) z = max{c T x : x e P C\ Z"}, 
where P C R™ is a given polyhedron, one often merely wishes to compute bounds 

z < z* < z 

on the optimal objective value z*. 

Apart from helping in implicit enumeration, upper and lower bounds can also 
provide important stopping criteria: 

• Having found a feasible point x with objective value c T x € \z_,z], one may 
decide to stop the algorithm if the gap z — z is smaller than a chosen error 
tolerance e. 

• In the extreme situation where e = 0, the bounds z,z provide a certificate of 
optimality. 

• Instead of predetermining an error tolerance s, one may bound the maximum 
computation time, and in this case one is interested in the gap z — z to esti- 
mate the quality of the best found objective value c T x as an approximation 
of the unknown value z*. 

Definition 6.1. A primal bound for problem (IP) is a lower bound z < z*. A 
dual bound is an upper bound z > z* . 

We remark that Definition 6.1 refers to the maximisation problem (IP). For min- 
imisation problems upper bounds are primal and lower bounds are dual. This con- 
vention guarantees that if a maximisation problem is reformulated as a minimisation 
problem, the notions of primal and dual bounds do not change. 

6.1. Finding Primal Bounds. The mechanism for finding primal bounds is 
conceptually simple: Any feasible point x G P fl Z" provides a primal bound c T x on 
the optimal objective value z*. 

Example 19. The IP 

z = max x\ — x 2 

s.t. x\ + 3x 2 < 5, 
xi, x-i > 0, 

x e z 2 

has x = (2,1) among its feasible solutions, hence 1 is a lower bound on z* (which 
takes the value 5 in this example). 

This may seem easy, but finding feasible solutions is often a difficult problem in 
itself. Typically, heuristics are used to overcome this problem. 
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EXAMPLE 20. A feasible solution to the 0-1 Knapsack problem 

max5xi + 8x 2 + 17^3 
s.t. 4xi + 3x 2 + 7x 3 < 9 

x e B 3 

can be found using the greedy heuristic, which is to grab objects that fit into the residual 
volume in order of decreasing value per unit volume. 

1. We note that €2/02 > c 3 /a3 > c\/a\. 

2. The residual volume is 9, so we grab object 2 and set x 2 = 1. 

3. The residual volume (after packing object 2 is 6 = 9 — a 2 . Object 3 does not 
fit, and we set x 3 = 0. 

4- The residual volume is 6, object 1 fits and we set x\ = 1. 



6.2. Finding Dual Bounds. There are two broad classes of methods to gener- 
ate dual bounds: 

i) Use a relaxation technique, such as LP relaxation, combinatorial relaxation, 
Lagrangian relaxation, semidcfinitc programming relaxation etc. (see the 
later parts of this course) . 

ii) Exploit duality. 

Definition 6.2. We say that the problem 

(RP) z R = max{/(x) : x e T C K"} 
is a relaxation of the integer programming problem 

(IP) z* = max{c T x : x e X C Z™} 

if the following two conditions are satisfied, 

i) X C T, 

ii) f(x) > c T x for all x G X . 



We note that these conditions imply 

z R = max fix) > max fix) > maxc T x = z* . 

Thus, if (RP) is an easier problem than (IP), then solving (RP) can be a useful 
technique for finding a dual bound for (IP). 

6.2.1. LP relaxation. An example of a relaxation technique we have already 
encountered is the following: 
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Definition 6.3. The LP relaxation of the IP 

(IP) max{c T x : x G PC)Z n } 
with formulation P = {x E K" : Ax < b} is the linear programming instance 

(LP) max{c T x : x G P}. 

The following trivial result shows that tighter formulations are indeed better than 
less tight ones, a result that we intuitively understood in Section 5 and that is now 
made precise: 

PROPOSITION 6.4. Let P 1 C P 2 C R n be two different formulations of the IP 
problem 

(IP) max{c T x : x G X}, 

that is, X = P 1 n IP = P 2 n IP, and let 

(LP1) z{ p = max{c T x : x G Pi}, 
(LP2) z% p = max{c T x : x G P 2 } 

be the optimal objective values of the LP relaxations corresponding to these two for- 
mulations. Then z\ F < z% F ', that is, the tighter formulation Pi yields a tighter dual 
bound. 

Proof. (LP2) is a relaxation of (LP1). □ 

We also say that (LP1) is a tighter relaxation than (LP2). 

6.2.2. Certificates of optimality. If a relaxation is tight enough it can either 
provide a certificate of infeasibility or optimality for the original problem, a situation 
we hope to encounter early on in a branch-and-bound approach: 

Proposition 6.5. Let 

(RP) z R = max{/(x) : x G T C M™} 

be a relaxation of the integer programming problem 

(IP) z* = max{c T x : x £ X C IP}. 

i) If x* is an optimal solution of (RP) for which x* G X and f(x*) — c T x*, 

then x* is an optimal solution for (IP), 
ii) If (RP) is infeasible, then (IP) is infeasible. 

Proof. 

i) Since x* G X, c T x* is a primal bound, and since x* is optimal for (RP), f(x*) 
is a dual bound for (IP). Therefore, c T x* < z* < f(x*) = c T x* shows that 
x* is optimal. 
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ii) Since X C T, we have T = => X = 0. 

□ 



Example 21. Consider the binary programming problem 

(BP) max 7xi + 4x 2 + 5x 3 + 2x 4 

s.t. 3xi + 3x 2 + 4x 3 + 2.T4 < 6 

x e B 4 . 

T/ie LP relaxation 

(LP) max 7a;i + 4x 2 + 5x 3 + 2x^ 

16I 4 

s.t. 3xi + 3x2 + 4x3 + 2.T4 < 6 
Xi <l (iG[l,4]) 
Xi>0 (ie[l,4]) 

/ms optimal solution x* = (1,1,0,0). Since x* is binary and we did not change the 
objective function, x* is also optimal for problem (BP). 



6.2.3. Combinatorial relaxation. Sometimes hard integer programming prob- 
lems can be relaxed by well-solved combinatorial problems such as those we study in 
later lectures. This is called combinatorial relaxation. 



Example 22. We will see that knapsack problems 

n 

(K) max ^ axi 

i=l 
n 

s.t. OiXi < b, 
x > 0, x e Z n , 

are often well-solved using dynamic programming in the case where the volumina di 
(i = 1, . . . , n) and b are integers. In the case where at and b are rational, (K) can be 
relaxed by the problem 

n 

(RK) max ^ gxj 

n 

s.t. ^\a,i\xi < [b\, 

i>0,iez", 

which falls into the well-solved category. 
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Example 23. Let (V,A) be a directed graph with arc weights Cij, (ij G A), and 
consider the travelling salesman problem 

n n 

min E 

CijXij 

i=l j=l 

S.t. ^ X i] = 1 (* G I 1 )"])) 

Y^ Xij = l (je[l,n]), 
X G B" x ". 

If we leave out the cut-set constraints the feasible set can only become larger, and we 
arrive at the assignment problem 

n n 

»=i j=i 

s.t. i I3 = 1 (i G [l,n]), 
^^• = 1 (je[l,n]), 
x G B" x ™. 

in /aier lectures we will see how the assignment problem can be solved efficiently. 

6.2.4. IP duality. A conceptually different approach to obtaining dual bounds 
is via a dual problem, if one can be identified: 

Definition 6.6. Consider an integer programming problem 
(IP) z* = max{c T x : x G X C Z"}, 

and let 

(D) s* = min{w;(ii) : u e U C R m } 

be a given second optimisation problem satisfying the relation 

c T x < w(u), Vx G X, u G U. 

Then (D) is called a weak dual of (IP). 

If furthermore z* = s* , (D) is called a strong dual of (IP). 

6.2.5. Further certificates of optimality. Sometimes dual problems can pro- 
vide certificates of infeasibility or optimality for the original problem: 
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Fig. 6.1. Example of an undirected graph with edge and node weights. 



Proposition 6.7. Let (D) be a weak dual of (IP), 
i) If x* G X and u* £ U are such that c T x* = w(u*), then x* is optimal for 
(IP)- 

ii) If (D) is unbounded (i.e., it has "optimal" value —oo), then (IP) is infeasible. 
Proof. 

i) c T x* is a primal bound and w(u*) a dual bound on z* . Therefore, c T x* < 
z* < w(u*) — c T x* shows that x* is optimal. 

ii) Since (D) is unbounded, for all Mel there exists u € U such that w{u) < M. 
However, if there did exist a x £ X, then there could be noneP such that 
w(u) < c T x. 

□ 

6.2.6. Examples of weak dual pairs. We begin by a trivial example which is 
nonetheless important, as it applies to all integer programmes: 



Example 24. The IP 

(IP) max{c T i : Ax < b, x E Z" } 

and the LP 

(D) mm{b T y : A T y > c, y e M™} 

form a weak dual pair. This follows trivially from the fact that (D) is the dual of the 
LP relaxation of (IP). 



The following is a less trivial example: 



Example 25. A matching in an undirected graph G = (V,E), where V is the 
set of nodes and E the set of edges, is a subset M C E of disjoint edges in E in the 
sense that no two edges share a common endpoint. For example, Figure 6.2 shows a 
matching of the graph of Figure 6.1. A covering of G is a subset R C V of nodes 
such that every edge in E has at least one endpoint in R. For example, Figure 6.3 
shows a covering of the graph from Figure 6.1. Now consider the problem of finding 
a maximum cardinality matching 

(M) max{|Af| : M is a matching} 
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Fig. 6.2. The red edges form a matching. 




Fig. 6.3. The red nodes form a covering. 



and the problem of finding a minimum cardinality covering 



min{|i?| : R is a 

RCV 




We claim that (M) and (C) form a weak dual pair. 

To prove this, let n := \V\ and m := \E\, and consider the n x m node-edge 
incidence matrix A = (aj ie ) which is defined as follows, 



For example, The node-edge incidence matrix of the graph from Figure 6.1 is given by 



Using the concept of node-edge incidence matrix, (M) can now be reformulated as 
follows, 



where 1 = [l ... l]. To see this, let Ai be the i-th row of A. Each component of the 
vector Xi corresponds to an edge in the graph. Now AiX < 1 and x 6 Z" implies that 
at most one of the Xi corresponding to an edge incident to node i is nonzero, and that 
this nonzero entry is then equal to 1 . Thus, all feasible x are binary and correspond 
to flag variables indicating which edges are selected, and these edges correspond to a 
matching. 




if node j is incident to (an endpoint of) edge e, 
otherwise. 



A = 



" 1 1 " 

1 1 

1 1 

1111 

1 

1 1 

.0 1.. 



(M) z* = max{l T x : Ax < 1, x € Z™}, 
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Similarly, (C) can be reformulated as 

(C) w* = mm{\ T y : A T y > 1, y e Z™ }, 

The argument to see this is similar, but here each component of y corresponds to 
a node in the graph. The LP relaxations of (M) and (C) are respectively the LP 
instances 

(P) z LP = max{l T a; : Ax < 1, x € M+ }, 
(D) w LP = min{l T ? y : A T y >l,y€ R™ }, 

which are duals of each other. By LP duality and the properties of relaxations, there- 
fore, for all x and y feasible for (M) and (C), 

l T x <z*< z LP = w LP < w* < l T y. 

This shows our claim. 

7. Total Unimodularity. Consider an integer programming problem 

(IP) max c T x 

s.t. Ax <b 
x>0 

x e Z", 

where we assume that A e Z mxn and b £ Z m have integer coefficients. It is natural 
to ask for whether we can recognise problem parameters (A, b) for which the LP 
relaxation 

(LP) max c T x 

s.t. Ax < b 
x>0 

actually solves the integer programming problem (IP) to optimality, rather than just 
providing a dual bound. Obviously, for this to happen there must exist an optimal 
solution of (LP) that is feasible for (IP), see Figure 7.1. 

A complete characterisation of such input data is not available, but we can give a 
partial answer in the case where the LP formulation is not degenerate, that is, where 
it has a unique optimal solution. In this case a basic feasible optimal solution must 
do the trick. In Section 2.2 we saw that such a solution is of the form 

%B = -^B ^' ^ N = 

for some nonsingular matrix As composed of the columns of [ A i ] that correspond 
to the indices in B, that is, As = CP for some submatrix C of [ai] and some 
permutation matrix P 3 . A first insight is the following: 

3 For the benefit of completeness, let us briefly recall the concept of permutation matrices. Let 
a : {!,...,«} — > {!,...,«} be a permutation, and let P = [pij] be the n X n matrix defined by 
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Fig. 7.1. The optimal solution of the IP on the left can be found by solving its LP relaxation. 
The IP on the right has a LP relaxation with infinitely many optimal solutions. One of these is the 
optimal solutions of the IP, but this solution is not found directly by the simplex algorithm, as basic 
feasible solutions correspond to vertices of the feasible polyhedron. 



Lemma 7.1. // | det(C) = 1, then (LP) solves (IP). 

Proof. We note that since det(P) = ±1, |det(vl B )| = |det(C)| • |det(P)| = 1. 
It suffices to prove that x* B = A^b is an integer vector, and since b has integer 
components, it suffices to show that Ag 1 = \J3ij] has integer components /3y. Let A % ^ 
be the (to — 1) x (to — 1) matrix obtained by deleting row i and column j from As- 
By Cramer's rule, 

_ det(Ag) Z 

□ 

A particularly lucky situation occurs when |det(C)| = 1 for all submatrices of 
[ A i ] that correspond to a set of basic variables. The following is a stonger notion 
that guarantees that this is the case: 

Definition 7.2. A matrix M is totally unimodular (TV) if every square sub- 
matrix (of any size) has determinant +1, —1 or 0. 

An immediate consequence of Definition 7.2 is that all components of a TU matrix 
M must be € {0, 1, —1}. To verify this, it suffices to consider submatrices of size lxl. 

Example 26. It can be verified by tedious evaluation of the determinants of all 



Pij = 1 if j = cr(i), and pij = otherwise. P is called the permutation matrix associated with a. 
If M is a m X n matrix with columns M- \, ■ ■ ■ , then M ■ P is the m X n matrix with columns 

M. a—in), • ■ ■ , M. a -ii n y In other words, right multiplication with P moves column i into the a(i)- 
th position. It is immediate to see that P T does the opposite, and hence, P T = P~ 1 . Therefore, 
det(P) 2 = det(P T ) ■ det(P) = det(I) = 1. By taking transposes it is trivial to see that ifMisanxm 
matrix with rows Mi,; , . . . , M n ,-. then P T ■ M is the nxm matrix with rows M^-i^ M CT -i ( n j . , 

i.e., row i moves into the <r(i)-th place. 
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square submatrices that the following matrix is TU, 

"0 1 0" 

1111 
10 111 
10 10 

1 

We understand that total unimodularity merely gives a sufficient criterion for 
when the integer programme 

(IP) max{c T x : Ax < b, x G Z™ } 
is solved by its LP relaxation 

(LP) max{c T x : Ax < b, x > 0}. 
The following result explains how restrictive this criterion is: 

Theorem 7.3. 

i) If A is TU, then (LP) solves (IP) for all b G Z m for which an optimal solution 
exists. 

ii) If (LP) solves (IP) for all b G Z m for which an optimal solution exists, then 
A is TU. 

Proof. Part i) follows from Lemma 7.1. Part ii) is left as an exercise. □ 

Corollary 7.4. Let A be TU and b G Z™. 

i) Let (D) be the dual of the LP relaxation of (IP). Then (D) is a strong dual 
of (LP). 

ii) The feasible set P = {x G W 1 : Ax < b, x > 0} of the LP relaxation is an 
ideal formulation of (IP), see Figure 7.2. 

Proof. 

i) Let w*, z LP and z* be the optimal objective values of (D), (LP) and (IP) re- 
spectively. From Lecture 2 we know that w* = z LP > z*, and from Theorem 
7.3 we know z LP = z* . Thus, w* = z*. 

ii) Every vertex of P can be turned into an optimal basic feasible solution of (LP) 
for an appropriate choice of the objective vector c. By Theorem 7.3 each of 
these vertices lies in the feasible set X = P n Z™ of (IP). By Theorem 4 hi) 
of Lecture 6, P is the convex hull of its extreme points, thus P C conv(X). 
Since on the other hand X C P, we have conv(X) C conv(P) = P. Thus, 
P = conv(X), as claimed. 

□ 

7.1. Turning Total Unimodularity into a Practical Tool. By now we re- 
alise that total unimodularity is a powerful concept, but verifying that a given matrix 
is TU seems a task of complexity exponential in the size of the matrix. To make this 
a useful practical tool we need simple criteria that help us recognise TU matrices. 
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Fig. 7.2. Totally unimodular constraint matrices yield ideal formulations of the underlying IPs: 
The feasible set of the LP relaxation is given by the convex hull of the IP-feasible points. 



The first type of useful criteria are rules by which small TU matrices can be as- 
sembled into larger ones. By applying the inverse of these rules, we may be able to 
recognise how to decompose a matrix into smaller parts whose total unimodularity is 
computationally cheaper to verify. The following result serve this purpose: 

Theorem 7.5. A m x n matrix A is TU if and only if any of the following 
matrices are TU, 

i) A T , 

ii) [A -A], 

Hi) A ■ P, where P is a n x n permutation matrix 

iv) P ■ A, where P is a m x m permutation matrix, 

v) [ j 2 J q], with Ji — ftfoo]^*' an d where I are identity matrices, blocks of 
zeros, and Pi,Qi permutation matrices of appropriate size. 

Proof. The proof is tedious but elementary: 

i) Let R C {1, . . . , m} and C C {1, . . . , n} be index sets, and let Arc be the 
submatrix of A obtained by deleting all rows not in R and all columns not 
in C. By definition, A is TU if and only if det(A RC ) e {0, 1, -1} for all R, C 
that correspond to square submatrices. Since 

det (Al R ) = det ((A RC ) T ) = det(A RC ), 

A is TU o A T is TU. 

iii) Let P be the permutation matrix corresponding to the permutation a : 
{1, .. .,n} -)■ {1, .. . ,n}, and let C = {jx, . . . ,j k }. Then 

(AP) R c = [A Rta -i {jl) ... A R ^-i Uk) ] . 

Let H — (hi, . . . , hk) be the ordered set {er -1 ^) : I = 1, . . . , k}, and let Q 
be the k x k permutation matrix such that 

[Ar.ct-iO'O ••• A B^-H Jk )} = A RH Q. 

Then det((AP) RC ) = det(A RH ) ■ det(Q) £ (±1) • det{A RH ). It follows that 
every square submatrix of AP has determinant in {0, 1, —1} if and only if 
every square submatrix of A satisfies this criterion. 
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ii) Every submatrix of A is a submatrix of [a -a]. Thus, if the latter is TU, 
then A is TU. For the converse, note that 

[A -A] RC =[A RCl -A RC2 ] 

for some index sets C\, C2 C {1, . . . , n}. If C\ n C 2 ^ 0, then the columns 
of [a -a] rc are linearly dependent, and the determinant is zero. Otherwise 
there exists a permutation matrix P, an index set C and a diagonal matrix 
D with diagonal entries ±1 such that 

l A - A ] RC = A RcPD, 

and then 



det( [A -A] ) = det(A Rd ) det(P) det(D) e {0, 1, -1}, 



since A was assumed to be TU. 

iv) This follows from a combination of iii) and i). 

v) Writing B = P? AQj , we find 



where 



'A Ji 




'Pi 


0" 




B I{ 




Qi 0" 


J 2 o_ 







P 2_ 




h 0_ 




Qi 





I 


0" 




h = 










'Pi 0" 




Q2 


" 


P 2 


5 








Since 



are permutation matrices, it follows from iii) and iv) that we can directly 
work with the matrix 



B h 
h 



Furthermore, if we can show the claim for matrices of the form 
only, then the claim follows from 



B TU &[B h] TU 



11 



TU 



B h 
h 



TU. 



B being a submatrix of [ s / ] , the total unimodularity of the latter clearly 
implies that the former is TU. If B is TU, then [b i]rc is of the form 



[BR,d Ji] 7 
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where Ji is of the same form as Jj, but possibly smaller. Since £r,Ci is 
TU, it thus suffices to assume that [b i] is square, show that det([s /]) € 
{0, 1, —1}, and apply the argument recursively. But if [b i] is square, then 
det([s /]) 7^ only if it is of the form 

~B 1 I" 

B 2 Oj ' 

and then 

dct([s /]) = ±dct(S 2 ). 

_B 2 being a square submatrix of B , and i? being TU , it follows that det ( [ b i ] ) 6 
{0, 1, —1}, as required. 

□ 



Let us now sec Theorem 7.5 in action: 



Example 27. The following matrix is TU, 

^10 -1 
11-1-1 

1 



Indeed, it is trivial to check that the matrix 
we find that 



1 
1 1 



is TU. By application of rule ii), 



10- 10 

11- 1-1 

is TU, and by application of v) it follows that 

10- 1 0' 

11- 1-1 
1 

is TU as well. By permuting the last two columns we find that A is TU. 

The second set of useful tools consist of sufficient criteria that can easily be 
checked and that allow us to identify some important families of TU matrices. The 
following theorem is a result of this kind, and further examples appear in the exercises. 

Theorem 7.6. Let A = [a^] be a matrix such that 
i) dij G {+1, — 1, 0} for all i, j. 

ii) Each column contains at most two nonzero coefficients, 



aij\<2 (je[l,n]). 



64 



Fig. 7.3. A network flow problem with production rates bi at each node i, arc rate capacities 
ha and rate costs ca . 



Hi) The set M of rows can be partitioned into (M\,M2) such that each column j 
containing two nonzero coefficients satisfies 



Then A is totally unimodular. 

Proof. The proof is by contradiction, assuming that A is not TU. Let B be a 
smallest submatrix of A such that det(B) ^ {0, +1, — 1}. Then all columns of B 
contain exactly two nonzero coefficients, for else there exists a permutation matrix P 
such that 



and then det(C) = ± dct(_B) ^ {0, +1, —1}, and C is the row permutation of a strict 
submatrix of B, contradicting the choice of B. Because of iii), adding the rows of 
B with indices in Mi and subtracting the rows with indices in Mi yields the zero 
vector, showing that the rows of B are linearly dependent and det(£>) = 0, which 
again contradicts the choice of B. □ 

7.2. Applications to Network Flows. The minimum cost network flow prob- 
lem is defined as follows: Given is a directed graph D = (V, A) with arc capacities hij 
for all G A. At each node i G V, bi units of liquid are produced per unit time 
(negative production is to be interpreted as consumption). The total consumption 
equals the total production. Liquid can flow via the network at rates up to the spec- 
ified capacities hij, but no spillage can occur. Pumping one unit of liquid per unit of 
time along arc G A incurs a costs cy. See Figure 7.3 for an example of such a 
network. The problem is now to determine which feasible flow will minimise the total 
cost. 

We will first formulate the min-cost network flow problem as an LP. Let Xij be 
the rate at which liquid is flowing along arc Let V + (i) = {j : G A} denote 
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i€M 2 



the set of nodes to which there is an arc from i (successor nodes), and V~(i) = {j : 
G A} the set of nodes from which there is an arc to i (predecessor nodes). We 
thus obtain the following LP model, 

min ^2 

CijXij 

(i,j)eA 

S.t. ^ x ij - J! X 3 i =bi ^ 6 V ) 

jev+(i) jev-(i) 
OKxij^hij ((i,j)eA). 



Before proceeding further, let us consider a specific example that will reveal the 
general structure of the constraint matrices of network flow problems: 



Example 28. The digraph (directed graph) with 3 nodes and arcs 
A = {(1, 2), (1, 3), (2, 1), (2, 3), (3, 2)} with capacities 3,2,4,3,6, transportation costs 
1, 1, 2, 3, 4 and node production rates 2, 4, —6 is modelled as 



min c T x 





' C ' 




' b ' 


s.t. 


-C 


x < 


-b 




I 




h 



x>0, 



where 



C = 



i i 

-1 o 
o -i 



-10 

1 1 -1 

0-11 



a=l2 
Xl3 

X21 
X23 
X31 
X32 



b = 



h = 



It turns out that the constraint matrix is TU. 

The findings of Example 28 remain true for general network flow problems: 



Proposition 7.7. The constraint matrix arising in a minimum cost network 
flow problem is totally unimodular. 

Proof. In complete analogy to iExample 28, the constraint matrix has the form 

" C ' 

-C . 
I 

It suffices therefore to show that C is TU. Each column of C has exactly two nonzero 
entries, one +1, the other one -1, because it enters in one constraint as an outgoing 
arc, and in another as an icoming arc. Therefore, the sufficient criterion of Proposi- 
tion I is satisfied with M 1 = M and M 2 = □ 

Corollary 7.8. In a minimum cost network flow problem, if the production rate 
vector b and the capacity vector h are integral, then 
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Fig. 7.4. The network of Example 28 



i) Each extreme point of the feasible polyhedron is integral (has integer coordi- 
nates). 

ii) If there exists an optimal flow, then there exists an integral optimal flow (since 

there exists an optimal extreme point), 
in) The constraints of the problem describe the convex hull of the set of integral 
feasible flows. 

The min-cost network flow problem is an LP, so the total unimodularity of its 
constraint matrix seems of no use at first sight. It turns out however that due to 
the unimodularity network flow problems can be solved via specialised algorithms 
(Ford-Fulkerson Algorithm) that are more efficient than the standard simplex method. 
Further, the advantage becomes immediately visible if we consider the shortest path 
problem, which is defined as follows: Given is a directed graph D = (V, A) with non- 
negative arc lengths for all (i,j) € A. Two nodes s,t € V are marked. The 
problem is to find a shortest path from s to t in D. 

Setting b s — 1, b± — — 1 and hi = for the remaining i s V, and 




1 if (i, j) is in the shortest path, 
otherwise, 



and hij = I for all arcs, we can model the shortest path problem as a special case of 
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Fig. 7.5. A shortest path problem interpreted as network flow problem, and an optimal solution 
revealed in bold arcs. 



the min-cost network flow problem with additional integrality constraints: 
(SP) z = mm ^2 

CijXij 

S.t. X s j — Xj s = 1 

jev+(s) jev-(s) 

j£V+(t) j£V-(t) 

jev+(i) jev-(i) 
0< Xij <l {(i,j)€A) 
x e l) A \. 

Total unimodularity now reveals the following useful properties of the shortest 
path problem: 

• Because of Corollary 7.8, (SP) is solved by a basic feasible optimal solution 
of its LP relaxation. 

• Due of Corollary 7.4, the dual 

(D) w LP — max 7r t — tt s 

S.t. TTj - 7Tj < Cij € A). 

of the LP relaxation is a strong dual of (SP). Setting tt s — without loss of 
generality, there is a very simple algorithm to solve (D) in &(\V\ x \A\) time 
(this is left as an exercise). 

8. Matroids and Submodularity. So far we know only one type of easy integer 
programming problem, those whose LP relaxation solves the problem. We will now 
discuss another type of easy problem, those that can be efficiently solved via greedy 
algorithms. 
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Fig. 8.1. Examples of a forest, a spanning forest, a tree, and a spanning tree. 



We start with the example of finding maximum weight trees and forests in undi- 
rected graphs, and then we generalise our findings to submodular optimisation prob- 
lems and the associated notion of matroids. 

Definition 8.1. Let G — (V,E) be an undirected graph with node set V and 
edge set E. 

i) A forest in G is a subgraph G' — (V',E f ) that contains no cycles, that is, any 

pair of nodes is connected in G' via at most one path in G' . 
ii) A tree in G is a connected forest, that is, a forest that contains a unique path 

between any pair of nodes (i,j). 
Hi) A spanning forest in G is a tree that covers every node in V. 

iv) A spanning tree in G is a spanning forest that is a tree. 

Proposition 8.2. Let G — (V,E) be an undirected graph with n = |V| nodes, 
the following characterisations are equivalent, 

i) G is a spanning tree, 

ii) G is a spanning forest that contains exactly n — 1 edges, 

Hi) G is an edge-minimal connected graph that spans V , that is, such that every 

node in V is covered, 
iv) G contains a unique path between every pair of nodes of V , 

v) G is a graph such that the addition of an edge not in E creates a unique cycle. 

Proof. Left as an exercise. □ 
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8.1. The Maximum Weight Spanning Tree Problem. Let G = (V,E) be 
an undirected graph with weights c e (e G E). The Maximum Weight Spanning Tree 
Problem is to find a spanning tree in G for which the sum of edge weights is maximised: 



(MWST) max l^c e :E'CE s.t. G' = (V, E') is a spanning tree 

leeE' 

Note that if G is not connected, then the feasible set is empty, and (MWST) is an 
infcasiblc problem. So far, (MWST) is a formulation as a combinatorial optimisation 
problem, but at the end of this section we will see that it also has a formulation as 
an integer programming problem that makes it a special case of a general family of 
IPs that can be solved via a greedy algorithm. Specialised to our problem, the greedy 
algorithm becomes the following: 

Algorithm 8.3 (Greedy Algorithm for (MWST)). 

1. Order the edges in E by nonincreasing weight c\ > C2 > • • • > c m , where Ct 
is the weight of edge et, set T° = 0. 

2. For t = 1 , . . . , m repeat 

i) If T l ~ x U {e t } contains no cycle, set T* = T t_1 U {e t }. Otherwise, 

it) If |T*| =n-l, stop with (V,T*) optimal. 

3. Stop with "no feasible solution". 



Example 29. Let G = (V, E) be the graph with node set V = {1, ... ,7}, edge set 

E = {{1,2}, {1,3}, {1,5}, {2,3}, {2,4}, {2, 5}, {2,6}, 

{3, 4}, {3, 6}, {4, 5}, {4, 7}, {5, 6}, {5, 7}, {6, 7}} , 

and edge weights c\.2 — 5, — 9, c^ 5 = 10 7 c 2 ,3 = 3 7 c 2 a = 3, c 2 ,5 = 5, c 2 ,6 = 6, 
c 3 ,4 = 12, c 3 , 6 = 9, c 4j5 = 2, c 4 , 7 = 4, c 5 , 6 = 9, c 5j7 = 6, c 6j7 = 8. The greedy 
algorithm chooses the following edges in order: {3,4}, {1,5} ; {1,3} ; {3, 6} ; {6,7}, 
{2,6}, see Figures 8.2-8.7. 

8.2. The Maximum Weight Forest Problem. This is a slight variation of 
(MWST) in which the aim is to find a forest in G = (V, E) of maximum weight, 

(MWF) max \^c e :E'<ZE s.t. G' = (V, E') is a forest i . 

le£E' J 

Note that, contrary to (MWST), (MWF) is always feasible. This problem is also a 
special case of the problem class we will discuss later, and hence it is also efficiently 
solvable via the greedy algorithm. Specialised to this case it takes the following form: 
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Fig. 8.2. Step 1 of the greedy algorithm applied to Example 29. 




Fig. 8.3. Step 2 of the greedy algorithm applied to Example 29. 



Algorithm 8.4 (Greedy Algorithm for (MWF)). 

1. Order the edges in E by nonincreasing weight c\ > C2 > • • • > c m , where c\ 
is the weight of edge et, and set T° = 0. 

2. For t = 1 , . . . , m repeat 

i) If c t < 0, stop with T 1 ^ 1 optimal. 

ii) If T l ~ l U {et} contains no cycle, set T* — T l ~ l U {e t }. Otherwise, 
Hi) If |T*| =rt — l or t = m, stop with T optimal. 
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Fig. 8.4. Step 3 of the greedy algorithm applied to Example 29. 




Fig. 8.5. Step 4 of the greedy algorithm applied to Example 29. 



8.3. Submodularity. Rather than proving that the greedy algorithm solves 
(MWST) and (MWF) efficiently, we derive this result for a much more general class 
of optimisation problems, see Theorem 8.12 below. 

Definition 8.5. Let n e N, N = {1, ...,n}, and let &>{N) := {A : A C N} be 
the power set of N , that is, the set of all its subsets. A function f : &(N) —> R is 
called 

i) submodular if 

f(A) + f(B) > f(A n B) + f(A U B) V A, B e &(N), 
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Fig. 8.6. Step 5 of the greedy algorithm applied to Example 29. Note that at this stage it is 
not edge {5, 6} that is chosen, despite its heavier weight in comparison to {6, 7}, because it would 
create a cycle. 




Fig. 8.7. Step 6 of the greedy algorithm applied to Example 29. After this step all nodes are 
covered and the algorithm stops. 



ii) and nondecreasing if 

f(A)<f{B) \/ACBe^(N). 



Definition 8.6. Let f be a nondecreasing submodular junction on &(N) con- 
structed so that /(0) = 0. The submodular polyhedron associated with f is defined 
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by 



P(f) = { 



x 



efix^O,^!^ f(S), VS g &(N) \ {0}}. 



Given a vector eel", the submodular optimisation problem associated with f and 
c is given by 



We note that (SOP) is an LP with 2" — 1 constraints. Therefore, even for mod- 
erately sized n the constraint matrix of (SOP) could not be held explicitly in the 
memory of a computer, and the simplex method could not be applied. It is therefore 
somewhat suprising that an efficient algorithm for solving (SOP) is nonetheless avail- 
able: 

Algorithm 8.7 (Greedy algorithm for (SOP)). 

1. Order the variables so that C\ > ■ ■ • > c r > > c r+ i > ■ • ■ > c n (where r = n 
ifc n >0), and let S" = 0. 

2. For i = l,...,r 

i) setS 1 = 

ii) setx i = f{S i )-f{S i - x ). 

3. For j > r, set Xj = 0. 



Theorem 8.8. The greedy algorithm solves the submodular optimisation problem 
to optimality in &(n(cj + logn)) time, where cj is the amount of work required per 
function evaluation of f. 

Proof. The complexity claim is immediate, as steps 2 and 3 take at most ncj 
time, and step 1 takes 0(n\ogn) time. Thus, the claim is all about the algorithm 
yielding an optimal vector x. 

Claim 1: x is (SOP)-feasible: Since / is nondecreasing, Xi = f(S l ) — /(S 14-1 ) > 
(i = 1, ...,r), and since Xi = (i = r + 1, . . . , n), the nonnegativity constraints 
x > are satisfied. Further, for all T e 0?{N) and j 6 T, applying submodularity to 
A = S^- 1 and B = 5^ ' n T, we find 



(SOP) max{c T x : x E P(f)}. 



f(S^ 1 ) + f{S^T)^f{A)+f{B) 

> f(A n B) + f(A U B) 

= f(S j - 1 nT) + f(S j ), (since j 6 T) 



and hence, 



f{Si) - fiSi- 1 ) < f(S* fl T) — fiSi- 1 n T). 



(8.1) 
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Thus, 

E x i = E (f( S ^ ~ /O^'" 1 )) ( b y construction of x) 

< Yl (/^nTj-J^nT)) (by inequality (8.1)) 

jeTr\S r 

< Y (f( sj n T) - / (S^ 1 n T)) (since / is nondccreasing) 
= f(S r nT)-f(Q))<f(T). 

Claim 2: The objective value at x, 

r 

»=i 

is optimal for (SOP): We prove this claim by constructing a dual feasible solution with 
the same dual objective value. The claim then follows from the LP duality theorem. 
It is easily verified that the LP dual of (SOP) is 

(D) min J2 f( S )-VS 

s-t. ys> Cj V j e N, 

y s >0 VSe&(N). 

For i = 1, . . . , r — 1 set j/s< = Cj — c i+ i. Set ysr = c r , and set ys = for all other 
S e &>(N). Then y is (D)-feasible, since y s > for all S e &>(N), and furthermore 
for j < r, 

r r—1 

E ys = ^2ys* = E( c * - c m) + c r = cj, 

S:jeS i=j i=j 



whereas for j > r, 



VJ y s = > Cj 

s-.jes 



The dual objective value at y is 

r r—1 

E = E /( si )( c ' - c *+i) + /(^k 

i=l i=l 
r 

= E C * (f(S i )-f(S i - 1 ))=z*. 



»=i 



We have seen that the submodular optimsation problem is an LP with 0(exp(n)) 
constraints that can be solved in 0(n). While this is in itself interesting, it is natu- 
ral to ask what the connections with integer programming are: Observe that if / is 
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integer valued, then the greedy solution of (SOP) is integral! We will now see that 
(SOP) is the LP relaxation of a related combinatorial optimisation problem. The 
above observation will then imply that this relaxation is tight. 



Definition 8.9. / : &{N) — > R is called a submodular rank function if the 
following are satisfied, 

i) f is submodular, 
ii) f is increasing, 
m) f (0) - 0, 

iv) for all S e &{N) and j <£ S, 

f(Su{j})-f(S)e{o,i}. 



Proposition 8.10. Let f be a submodular rank function. Then 

i) f(A) < \A\ for all A e &>{N), 
ii) Iff (A) = \A\ then f(B) = \B\ for all B C A. 
Hi) Let x A be the incidence vector of A £ 0?{N), defined by 




if J G A, 
otherwise. 



Then x A G P(f) if and only if f(A) = \A\. 

Proof. The only nontrivial part is hi). If f(A) < \A\ then J2jeA x f = 1^1 > f(A) 
and x A £ P(f). Conversely, if f(A) = \A\, then for all S £ &>(N), 

□ 



Definition 8.11. Let f : &>{N) :— > N be a submodular rank function, and let 
each element j of N be given weight Cj . 

i) A set A e 3^{N) is called independent under f if f(A) = \ A\. 

ii) The set & of independent subsets of N under f is called the matroid defined 

byf. 

Hi) The maximum weight independent set problem associated with f is 

(MLSP) max \ ^ ■ Ae& 
[jeA 



Theorem 8.12. The greedy algorithm solves (MLSP). 
Proof. By Theorem 8.8, the greedy algorithm solves 

(SOP) max i c i x i : x e P U) 
[ 3 eN 
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Since / is a rank function, the greedy solution x* of (SOP) is binary, x* e B™. 
Furthermore, 



(MISP) max 




where the last equivalence follows from Proposition 8.10 iii). This shows that (SOP) 
is the LP relaxation of (MISP). Since the optimal solution x* of (SOCP) is (MISP)- 
feasible, it follows by relaxation that x* is also (MlSP)-optimal. □ 

Example 30. Let G = (V,E) be an undirected graph, and let 

f : &{E) -> R, 

A h-> max{|F| : F C A, (V, F) is a forest}. 

We claim that f is a submodular rank function. The associated matroid & is the set 
of forests in G. Thus, the associated maximum weight independent set problem is the 
same as the maximum weight forest problem for G. To prove our claim, we first see 
that for any edge set A C E, f(A) = \ A\ if and only if A is a forest, so & is indeed 
the set of forests in G. We also clearly have /(0) = 0, and 



for all B C A (since forests in B are then also forest in A). Further, let A = B U {e} 
with e ^ B, and let F be a forest of maximum cardinality in A. Two cases can 
occur: either there exists a forest H C B such that \H\ — \F\, and then (8.2) implies 
f(B) = f(A), or else all cardinality maximising forests in A contain e. But then 
F \ {e} is a forest in B of cardinality \F\ — 1, which must be of maximum cardinality 
by nonexistence of H and (8.2). Thus, in both cases we have 



So far we have shown that f is an increasing rank function. It remains to show 
submodularity. This will be part of the problem sheets. 

Note finally that applying Algorithm 8.7 to the submodular optimisation problem 
associated with f is the same as Algorithm 8.4- 

9. Augmenting Path Algorithms. In this section we investigate matching 
problems which can be efficiently solved using augmenting paths. Recall the notions 
of matching and covering introduced in Example 25: 

Definition 9.1. Let G = (V,E) be a graph. A matching M C E is a set of 
disjoint edges, that is, such that no two elements of M share a common node as an 
endpoint (see Figure 6.2 for an example). A covering by nodes R C V is a set of 



f(B) < f{A) 



(8.2) 



f(A)-f(B)e {0,1}. 
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Fig. 9.1. Examples of an alternating and an augmenting path. 

nodes such that every edge in E has at least one endpoint in R (see Figure 6.3 for an 
example). 



In Example 25 we showed that the maximum cardinality matching problem 
(M) maxjlMl : M is a matching) 

MCE 

and the minimum cardinality covering problem 

(C) min{|i2| : R is a covering) 

RCV 

form a weak dual pair, that is, \M\ < \R\ for all matchings M and coverings R. Our 
goal is now to build an algorithm for a special case of the max cardinality matching 
problem that finds successively better matchings. The following concepts are used to 
find updates to a matching M of a graph G = (V, E): 



Definition 9.2. An alternating path with respect to M is a sequence of nodes 
P = (uq, v\, . . . , v p ) such that 

i) e i+ i := (vi,v i+ i) € E for (i = 0, . . . ,p - 1), 

ii) e$ 6 M for i even, and e$ 6 E \ M for i odd, 
Hi) Vq is exposed (i.e., not covered by any e € M), 

Furthermore, P is called augmenting path if p is odd andv p is exposed, see Figure 9.1 



Lemma 9.3. If P is an augmenting path with respect to M then 

M' = (M U P) \ (M n P) 

is a matching with \M'\ — \M\ + 1, where we write M U P for M U {(vj, : i = 

0, ...,p-l} and M n P for M l~l {(«*, v i+ i) : i = 0, . . . ,p - 1}. 

Proof. Since vo and v p are exposed in M, they are incident to exactly one edge in 
M'. Any other node along P remains incident to exactly one edge, although this edge 
changes from an even one to an odd one. For any node that doesn't lie on P any inci- 
dent edge lies in M iff it lie in M'. Thus, all nodes are incident to at most 1 edge in M', 
so that M' is a matching. Finally, MnP = {(uj, f j+i) : i = 1, 3, . . . ,p — 2} is of cardi- 
nality (p-l)/2 and M'nP = {(v i7 v i+ i) : i = 0, 2, . . . ,p- 1} of cardinality (p+l)/2. 
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Fig. 9.2. Two matchings (one red, one blue) of a graph and their difference graph. 



Since M and M' coincide outside of P, \M'\ = \M\ - (p- l)/2 + (p+l)/2 = \M\ + 1. □ 

The following result now provides an optimality condition for the maximum car- 
dinality matching problem: 

Theorem 9.4. A matching M for G — (V, E) is of maximum cardinality if and 
only if there does not exist an augmenting path with respect to M . 

Proof. The "only if" part follows from Lemma 9.3. To show the "if" part, we 
assume that M is not of maximum cardinality and construct an augmenting path. 
Let M' be a matching such that \M'\ > \M\, and consider the graph (V,F) where 
F = (MUM') \ (MnM 1 ). Since each node v G V is incident to at most one edge from 
M and M' each, it is incident to at most two edges in F. The connected components 
of (V,F) therefore consist of paths and cycles, see Figure 9.2. Since M and M' are 
matches, the edges in each connected component must alternate between elements of 
M and M' . For a component that is a cycle, this can only be true if it contains an 
equal numbers of elements from M and M' . Likewise, a component consisting of a 
path of even length contains equal numbers of elements from M and M' . Since F 
contains more elements from M' than from M, at least one of the connected com- 
ponents must be a path of odd length starting and finishing with an edge from M' . 
This path is an augmenting path. □ 



Theorem 9.4 leads to a conceptual algorithm for the max cardinality matching 
problem: 

79 




Algorithm 9.5. [Conceptual algorithm] 

1. Set M := 0. 

2. While 3P augmenting path w.r.t. M, repeat 

i) find P 

ii) M <- (M UP)\ (M n P) 



9.1. Maximum Cardinality Matching in Bipartite Graphs. To make this 
a practical algorithm, we need to be able to efficiently find an augmenting path when 
one exists, or else to find a certificate of inexistence. In bipartite graphs this is possible: 

Definition 9.6. A graph G = (V,E) is bipartite if the set of nodes can be split 
into two disjoint parts V = Vi U V% such that every edge e £ E is incident to a node 
from V\ and one from V2, see Figure 9.3. 

Since any augmenting path contains an odd number of edges, one of its end nodes 
must lie in V\ and the other in V2. To find an augmenting path w.r.t. a matching, it 
therefore suffices to start tracing it from V\. The following algorithm systematically 
builds up the beginnings of all possible augmenting paths until one is identified or a 
certificate of nonexistence is found. It employs a system of labelling (useful to trace 
the path back from its end to its starting point) and scanning (to keep account of the 
nodes that have already been looked at). Before reading the following algorithm, the 
reader is advised to study the example of Figure 9.4. 

Algorithm 9.7. [Finding an Augmenting Path in a Bipartite Graph] Given are 
a bipartite graph G — (Vi, V2, E) and a matching MCE as inputs. 

1. Initialisation: Label each exposed node in V\ with *. 

2. While there are unscanned labels, choose a labelled unscanned node i, and do 
the following. 
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Fig. 9.4. In pursuit of an augmenting path with respect to the matching given by the red edges, 
we first label all the exposed nodes in V\ (in this case only node 2). Since an augmenting path 
starting in node 2 would have to pass over to V2 via an unmarked edge, we label all unlabelled 
candidate nodes (in this case only node 5) by 2. To continue an augmenting path, one has to pass 
back to Vi via marked nodes (thus labelling node 3 by 5), and then again to V2 via an unmarked 
edge. If, carrying on this fashion, an exposed node in V2 is found, an augmenting path will have 
been traced out by the labels. 

i) if i G V\ and 3j 6 V% unlabelled s.t. € E \ M , label j with i, 
ii) if i G Vi and jBj G V2 unlabelled s.t. (i,j) G E\M , mark i as "scanned" , 
Hi) ij ' i G V2 exposed, construct an augmenting path P by tracking labels back 
to * and stop, 

iv) if i G V2 not exposed, find the unique edge (j, i) G M, label j with i and 
mark i as "scanned" . 
3. Stop with certificate of nonexistence of an augmenting path. 



Theorem 9.8 (Correctness of Algorithm 9.7). Algorithm 9.7 finds an augmenting 
path with respect to M or detects nonexistence thereof in &{m) time, where m = \E\. 

Proof. The bound on the running time arises from the fact that every time a node 
is scanned, a different edge is explored. Hence, step 2 takes are at most m iterations, 
each at a constant cost. It is trivial to see that when the algorithm enters step iii), an 
augmenting path P is found. It suffices thus to prove that when the algorithm enters 
step 3, then an augmenting path does not exist. By Theorem 9.4, this is the same 
as claiming that the input matching M is of maximum cardinality. It follows from 
Example 25 (weak duality between max cardinality matching and min cardinality 
covering problems) that to show this, it suffices to construct a covering R C V such 
that \M\ = \R\. 

Let Vy, V2 be the nodes of V\ and V% that are labelled by the time the algorithm 
reaches step 3, and let , V 2 ~ be the unlabelled nodes. We claim that R = V{~ U V 2 + 
is a node covering of E such that \M\ = \R\. The only way R could fail to be a node 
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covering is for there to exist an edge G E with i G and j G . There 

cannot exist an edge (i,j) G E\M with i G and j G KT, since all nodes have 
been scanned, and when i was scanned j would have been labelled. There also cannot 
exist any edge (i, j) G M with i G and j G V 2 ~ , since then i would not be exposed, 
and the only way it could have got labelled would be if j got scanned. Therefore, all 
edges in E are covered by nodes from R, as claimed. 

It remains to argue that \M\ = \R\: By assumption of termination in step 3, no 
exposed nodes have been found in V 2 + . For every j 6 V 2 + there therefore exists i G V\ 
with (i,j) G M. When j got scanned, i was labelled, so we have i G V± . On the 
other hand, for every i G V{~ there exists j G V2 such that (i,j) G M, for otherwise 
i is exposed and would have been labelled with *, so that i G V-^~ . It must further 
be the case that j G , for otherwise i would have been labelled with j when j was 
scanned. Every element in V{ U V 2 + can thus be associated with a different edge from 
M, so that we have 

\R\ = \V 1 -UV+\<\M\. 
But weak duality says that \M\ < \R\, so that we have \M\ = \R\, as claimed. □ 

Corollary 9.9 (Konig's Theorem). In bipartite graphs the maximum cardinality 
matching problem and the minimum cardinality covering problem are strong duals of 
each other. 

Proof. Given a bipartite graph, consider running Algorithm 9.7 repeatedly, start- 
ing with the empty matching M = and using any augmented paths found to find a 
better matching. The algorithm must terminate in finite time, as the cardinality of 
a matching can increase at most m times. In the last iteration, Algorithm 9.7 must 
terminate in step 3. The proof of Theorem 9.8 then shows that there exists a covering 
with the same cardinality as the matching that has been constructed. □ 

Definition 9.10. Any labelling of G that can be produced by a valid partial run 
of Algorithm 9.7 is called a valid labelling with for G. 

We remark that if a valid labelling is already available, then the initialisation step 
1 of Algorithm 9.7 can be skipped, and the scanning continued in step 2. We refer 
to this situation as a warm start. Using Algorithm 9.7 as a subroutine in step i) of 
Algorithm 9.5, we arrive at a polynomial-time algorithm for the bipartite maximum 
cardinality problem: 

Algorithm 9.11. [Max Card. Matching Algorithm for Bipartite Graphs] Given 
is a bipartite graph G = (Vi, V~2,E) as input. 

1. Initialisation: Choose an initial matching M (e.g., via a greedy heuristic or 
by setting M = %). 

2. Repeat 

i) Run Algorithm 9.7 on the input G, M . 

ii) If the output of Algorithm 9.7 is an augmenting path P, update 



M <r- (M UP) \ (MHP). 
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Fig. 9.5. The greedy inpection of the graph of Example 31 reveals a matching of cardinality 9. 

Hi) Otherwise, if the output of Algorithm 9.7 is "no augmenting path exists", 
stop with M optimal. 



Example 31. Ten researchers are engaged in a set of ten projects. Let Si denote 
the researchers working on project i for i = 1,...,10. To keep track of progress 
or problems, management wishes to designate one person working on each project 
to report at their weekly meeting. Ideally, no person should be asked to report on 
more than one project. Is this possible or not, when S\ = {3,7,8,10}, S2 = {4,8}, 
S 3 = {2,5,7}, Si = {1,2,7,9}, S 5 = {2,5,7}, S 6 = {1,4,5,7}, S 7 = {2,7}, S 8 = 
{1, 6, 7, 10}, S 9 = {2, 5}, Sia = {1, 2, 3, 6, 7, 8, 10} ? 

To answer this question, we set up a bipartite graph with nodes V\ = {1, . . . , 10} 
corresponding to researchers, and V2 = {1', ■ ■ ■ , 10'} corresponding to projects. The 
edges are defined as follows, 

(i,j r ) e E & i e Sj. 

A maximum cardinality matching of this graph then gives a maximum number of dif- 
ferent researchers to present projects. By inspection we find a matching of cardinality 
9 is available, see Figure 9.5. Applying Algorithm 9.11, we find that we run out of 
labels to scan before an exposed node in V2 is identified, see Figure 9.6. By Theorem 
9.8, this proves that a matching of cardinality 10 does not exist, and hence that there 
is no feasible assignment of speakers that would see a different person presenting each 
project. 

Theorem 9.12 (Correctness of Algorithm 9.11). Algorithm 9.11 solves the max- 
imum cardinality matching problem on a bipartite graph G — (Vi,V2, E) in &{m 2 ) 
time, where m = \E\. 

Proof. The correctness of Algorithm 9.11 is immediate from Theorems 9.4 and 
9.8. The bound on the running time arises as each iteration of step 2 adds a new edge 
to M, and hence there are at most m iterations. Each iteration necessitates a run of 
Algorithm 9.7, which takes &(m) time by Theorem 9.8. □ 
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Fig. 9.6. Labeling and scanning reveals that the greedy matching was optimal. 



9.2. The Hungarian Algorithm for the Assignment Problem. In Sec- 
tion 1 we had a brief encounter with the assignment problem 

n n 

(AP) max ^2 X! 

CijXij 

i=i i=i 

n 

s.t. }]xij = 1 (i e [l,n]), 
i=i 

n 

^xy=i (je[l,n]), 

i=l 

lij-eB (i,j€[l,n]). 

By now we recognise that this is an easy type of IP, since its constraint matrix is 
easily recognisable as satisfying the sufficient criterion of Theorem 7.6 and hence as 
a totally unimodular matrix. Thus, the assignment problem can be solved via its LP 
relaxation 

n n 

(P) max ^2 

CijXij 

1=1 j=l 
n 

s.t. Xj j = 1 (i £ 

3=1 

71 

^2x tj = l [je[l,n]), 

i=l 

%ij > («, j G 

We will now discuss a specialised algorithm based on the maximum cardinality match- 
ing problem that solves the assignment problem in polynomial time. This method is 
known as the Hungarian Algorithm. 
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We start by noting that the dual of (P) is 



n n 

(D) min u * + EI v 3 

i=i j=i 

s.t. Ui + vj > Cij = l,...,n). 

Since the assignment problem is neither infeasible nor unbounded, the same must 
be true for its LP dual (P), as well as for the dual (D). Thus, strong LP duality 
applies. Let us now derive a slightly weaker form of strong duality independently in 
the specific context of problems (P) and (D). The insight gained from this derivation 
will be useful in the construction of the algorithm. 

The following observation holds both for the assignment problem and the TSP 
problem: 



Lemma 9.13. Given real numbers {wi}™ =1 and {vj}" =1 , let us consider the ana- 
logue of problem (P) with changed weights Cij — — Ui — Vj, 



n n 



(P') max EI EI Xij 

i=l j=l 
n 

s.t. ^Xij = 1, (ie [l,n]), 

3=1 

n 

E^J = C? e t 1 ."]), 

i=l 

Xij > 0, e [l,n]). 

Then the objective value of (P') differs from that of (P) only by an additive constant. 
Proof. For all x feasible for (P) (or, equivalently, for (P')), we have 

n n n n 

EI EI Cij x ij = EI 5Z( C *J — U i ~~ V j) X ij 
i—l j=l z— 1 j — 1 



n n 



= EEw -J2 v i EI : 

»=i j=i i=i V j=i / i=i Vi=i 



n n 



= XI EI ^j^' - EI u * ~ EI ^ t 9 - 1 ) 

z— 1 j — 1 i— 1 j — 1 



where we have used Y^j=i x ij = 1 an d Sr=i x *j = 1- ^ 



The above observation leads to the following optimality criterion: 



Theorem 9.14. If there exist vectors u,v E R" and an (AP)-feasible assignment 
x* £ B" 2 such that 

i) c^ = c^ — Ui — Vj < for all and 
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ii) x*j = 1 => Cij = 0, 
then x is an optimal assignment and has objective value Y^ii=i Ui Sj=i v j- 
Proof. Condition i) implies that for any feasible assignment x we have 

Condition ii) implies that the assignment x* satisfies j CijXij — 0, and by (9.2) 
this implies that x* is optimal for the assignment problem with changed weights c. 
By Lemma 9.13 it follows that x* is also optimal for the assignment problem with 
original weights c, and that it has the claimed objective value. □ 

Remarks: 

i) Note that condition i) simply says that (u, v) is dual-feasible, whereas con- 
dition ii) says that strict complementarity holds. Theorem 2.8 showed that 
these two conditions together are necessary and sufficient conditions for LP 
optimality. Here we found a weaker result by only claiming sufficiency. 

ii) If the dual optimal values of (u*,v*) were known, the primal optimal x* would 
only take nonzero values on the "residual" graph (Vi,V2,E) with edge set 

E = {(i,j) : cy =0}, 

where Cjj = — u* — u*. As long as we could choose an x* with support 
on this graph and satisfying the assignment constraints, this vector would 
correspond to an optimal assignment. But choosing an x* that satisfies the 
assignment constraints is the same as finding a matching of size n, and this 
can be found using Algorithm 9.11. 

The latter remark suggests an algorithmic approach to solving the assigment prob- 
lem: Guess (u*, v*), compute the residual graph, and solve the maximum cardinality 
matching problem. If the optimal matching is of cardinality n, we found an optimal 
assignment. Otherwise, we need to improve our guess of (u*,v*). We call solving a 
max cardinality matching problem a primal step, since it concerns the computation 
of the primal variables x. Correspondingly, the updating procedure of (u, v) will be 
called a dual step. As we will see, the final labelling information from the primal step 
preceding a dual step will provide the information on how to update the dual variables. 

Algorithm 9.15. 

1. Initialisation: Find dual feasible (u,v), set M — 0, label all nodes in V\ with 
* (exposed), set = V\, V 2 + = 0, and take a primal step. 

2. Primal Step: 

i) Compute Cij = Cij — Ui — Vj, (i,j — l,...,n), find the edge set E = 

■ Cij — 0} and the subgraph G = (Vi,V~2,E). 
ii) Mark the nodes in as unscanned and warm-start Algorithm 9.11 
(bipartite maximum cardinality matching) with the inherited labelling, G 
and M as inputs. Remove the labels only when an augmenting path is 
found. Call the final matching M* . 
Hi) If \M*\ = n then M* corresponds to an optimal assignment; stop, 
iv) If \M*\ < n, replace M by M* , and go to 3. 
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3. Dual Step 

i) Compute 5 := min{|cjj| : i E V±~,j 6 V 2 ~}. 

ii) Replace m <— Ui — 5 (Vie V^), and Vj Vj + S (Vj G V" 2 + ), and go to 
2. 



Before we analyse this algorithm further, let us see it in action for a specific ex- 
ample. 

Example 32. Find a lower bound on the optimal value of a 6-city travelling 
salesman problem instance with distance matrix 
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by solving its assignment problem relaxation. 

We start by setting up the corresponding assignment problem graph. For this 
purpose, we split each node i of the TSP into two nodes i and i' and set up a complete 
bipartite graph as follows: Arc (i,f) in the new bipartite graph is given the edge weight 
of the arc from the TSP. Arcs which don't have a corresponding arc in the 

TSP and should be avoided are given a weight corresponding to an a priori upper 
bound on the length of a tour, e.g., 



Civ — nmaxcjj + 1 = 73, 
so that these arcs cannot appear in an optimal assignment solution. The TSP 



i=i j=i 
s.t. ^2 x ij = 1 (* e 

J2 x ij = 1 (je[i,n]), 

£5>«>1 (ScV,S^9), 

ieS j(£S 

x e B nxn . 
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can now be relaxed by the assignment problem 



i=i j>=i 

s.t. ^2 X H' = 1 (* G I 1 '"])' 

Xij> = i C?'e[i,n]), 
x e B" x ". 

Algorithm 9.15 is designed for a maximisation problem instead of a minimisation 
problem. Thus, we multiply the objective function by — 1 and find the vector x* that 
maximises the problem 



max Y 

Cij'Xiji 

»=i j'=i 
s.t. ^2 x ij' = 1 (i e 



i (ie[MD, 



wrei/i the weights matrix 
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Having completed our preparation works, we are now ready to start Algorithm 9.15. 

Initialisation: Find dual feasible variables (u,v). One possible solution would be 
to set (u,v) = (0,0). However, we can do better, as we want to introduce as many 
edges into the residual graph as possible: Setting Vj = max^ Cij, c~ij = Cij — vj, we find 



and 
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[ C ij] = 



Our choice of v guarantees that each column contains at least one zero. We can further 
guarantee that each row contains at least one zero too, by choosing Ui = maxj Cjj , 



Fig. 9.7. The first residual graph of Example 32. 



T 



-1 -8 

-2 

-9 -5 

-4 -3 ' 
73 

-2 -67 

The residual graph is now given in Figure 9. 7. 

Since the maximum cardinality matching algorithm builds up a matching by suc- 
cessive enlargements via augmenting paths, it is advanageous to start this procedure 
with a good initial matching. Such an initial matching can be chosen greedily, for ex- 
ample. By inspection we find the greedy matching M — {(1, 4'), (2, 5'), (5, 6'), (6, 2')}. 
Now using the labelling algorithm to find augmenting path, we obtain the lablelling 
of Figure 9.8. No exposed node is found in V2, hence the greedy matching was of 
maximum cardinality. No augmenting path exists! 

Since the matching found so far is of cardinality less than 6, we need to improve 
our guess of the dual variables by taking a dual step: First we find 

S = mm{-Cij : i G {1, 3, 4, 6}, j G {l', 3', 5', 6'}} = min 

And then we compute the new edge weights, 

c l3 ^c l3 (i G V+,j G V+), 
Ct 3 ^Ci 3 (i G V{~,j G Vjj - )) 
Cij-^Cij-5 (i G Vf,j G Vf), 
c i:i <- Cij +6 (ie Vf,j G Yf). 
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Cij = c^ — Ui, so that 



and 



u = 
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Fig. 9.8. The labelling of the first residual graph, proving that the greedy matching was of 
maximum cardinality. 




Fig. 9.9. The second residual graph. Note that none of the edges already in the matching have 
been removed, and that the old labelling remains valid, but that the labelled nodes in V\ now need 
rescanning, as they are incident to new edges that are not in the matching. 



Thus, the revised edge weight matrix is as follows, 
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Using this matrix, we can construct the revised residual graph of Figure 9.9. 

The old labelling remains valid and can be reused in the next primal step. Res- 
canning node 1 identifies an augmenting path 3 — ^ 4' — > 1 — >■ 3', see Figure 9.10. We 
update the matching and start labelling again, see Figure 9.11. This leads to a further 
augmenting path 4 — > 1', and to the optimal assignment of Figure 9.12. 

90 



Extracting the optimal decision vector x* , we find 

Xij = 1 & G {(1, 3), (3, 4), (4, 1)} U {(5, 6), (6, 2), (2, 5)}. 

This corresponds to the two disjoint subtours 1 — > 3 — > 4 — > 1 and 5 — > 6 — > 2 — > 5 
iwt/i ioia? length 2 + 3 + 4 + + 3 + 4= 16. Thus, 16 is a lower bound on the shortest 
TSP tour. 

Let us now revert to the theory and convince ourselves that Algorithm 9.15 pro- 
duces the correct output. To do this, we proceed via a series of Lemmas: 

Lemma 9.16. Whenever a primal step terminates with \M*\ < n, it is the case 
that |Vi + | > |K> + |. 

Proof. Since the last iteration of the matching algorithm failed to find an exposed 
node in V2, all labelled nodes j e are matched by (i,j) G M* with some i E 
which was labelled with j. Thus, 1^+ > In addition, since \M*\ < n there is 

at least one exposed node i e V + , labelled with *, so that |V^ + | > \V^\. □ 

Lemma 9.17. In every dual step we have 5 > 0. 

Proof. Lemma 9.16 shows that x ^ 0. Further, there do not exist i £ 
and j G V2 such that cy = 0, for else G -B and j would have been labelled 

when i was scanned. By dual feasibility of u,v, it follows that < for all i G V^, 

Lemma 9.18. 

5iep amounts to setting 

<- (i G Vi~,j G V^ + ), 
Cij <- (i G Vf, j G y 2 _ ), 
Cij <- Cij - (5 (i G Vf , j G V^), 
cij <- Cij + 5 (i e j eVf). 

ii) The new (u,v) obtained in each dual step is feasible for (D). 
Proof. By the choice of S, this is trivial to check. □ 

As a consequence of this lemma, the new residual graph contains at least one new 
edge with i G and j G V^ - , and since M C x , all the edges of M 
appear in the new residual graph. The process of finding a matching of size n can 
thus be continued. 

Keeping the labels of and V£ , and marking the nodes of as "unscanned" 
yields the labelling that would be obtained if the nodes of the dual-step-updated 
graph G new — (Vj.,V2,E new ) were scanned in the same order as in the old graph 
qoU = (y 1 ,V 2 ,E old ). Note that all the edges used under the previous labelling are 
now in x V^f and thus remain in the new graph. The nodes in need to be 
rescanned because they can be incident to new edges leading to nodes in V 2 ~ . 
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Lemma 9.19. 

i) In each primal step in which no augmenting path is detected the cardinality 
of V 2 + increases by at least one from the labelling information inherited from 
the previous primal step. 

ii) After at most n dual steps, \M\ will increase in the next primal step. 

Proof, i) Since is inherited from one primal step to the next, |V 2 + keeps 
increasing until an augmenting path is detected. The increase is strict because there 
exists E x V 2 ~ that will have been added to the residual graph in the dual 
step in-between, so that j joins upon rescanning i. 

ii) An augmenting path is detected at the latest when |V 2 + | > \M\. By part i), it 
takes at most \M\ + 1 primal steps before this happens. □ 

Theorem 9.20. Algorithm 9.15 solves the assignment problem in &(n 3 ) work. 

Proof. The final output of Algorithm 9.15 corresponds to an optimal assignment, 
because the vector defined by 

Xij = 1 e M 

together with the final (w, v) satisfies the criteria of Theorem 9.14. In fact, < by 
dual feasibility of (u, v), and x i3 = 1 => (i,j) e E c i3 — 0. Since \M*\ < n, \M*\ 
can increase at most n times. By Lemma 9.19, it takes at most n primal and dual 
steps between increases. Due to the re-use of the labelling information, the overall 
complexity of finding an augmenting path is 0{n 2 ). Thus, the overall complexity is 
&{n-n 2 ). □ 

We remark that if the labelling is computed from scratch every time the bipartite 
maximum cardinality matching algorithm is called in a primal step, the complexity 
is i^(n 4 ): Each such call has complexity &(n 2 ), and by Lemma 9.19 there are @(n) 
calls between increases in \M\, and it takes n increases until \M\ = n. 



10. Dynamic Programming. In this section we will discuss an approach by 
which some combinatorial optimisation problem can be solved by recursively reduc- 
ing the problem dimension. This is called dynamic programming. We illustrate this 
technique with the uncapacitated lot-sizing and integer knapsack problems, but the 
idea is much more general and can be applied in many contexts. 



10.1. Uncapacitated Lot-Sizing. The uncapacitated lot-sizing problem con- 
cerns the following situation: A factory produces a single type of product in periods 
1, . . . , n. Production in period t incurs a fixed cost f t , and each unit produced costs 
p t . Earlier production can be kept in storage for later sale at a unit storage cost of h t 
during period t. Initial storage is empty. Client demand in period t is d t . See Figure 
10.1 for an illustration. How should we plan production so as to satisfy all demands 
and minimise total costs? 
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Fig. 10.1. The uncapacitated lotsizing problem. 




Fig. 10.2. The lot-sizing problem as a network flow. 



In Example 7 we found the BP model 

n n n—1 

min y^p t x t + ^ fty* + ^ h t s t 
x,v,s t=i t=i t=i 

s.t. s t -i + x t — d t + s t (t=l,...,n), 

x t <My t (t=l,...,n), 

so = s n = 0, s 6 K™ +1 , x e M™, y e B™. 

where M is some a-priori bound on the total consumption 52™=i ^t- 

ULS can also be interpreted as a min-cost network-flow problem with fixed charge 
as follows: Consider the directed graph G — (V, E) with V = {0, 1, . . . , t, . . . , n} and 
E = {(0, t) : t = 1, . . . , n} U {(t,t + 1) : £ = 1, . . . , n - 1}. Node is a source with 
production D = Y^i=i^i- Node t has production — d t (t = l,...,n). See Figure 
10.2 for an illustration. Now let xt be the flow along (0, t) for t = 1, . . . , n at a cost 
/t + ptXt if xt > 0, and zero otherwise. Let further St be the flow along (t, t + 1) for 
t — 1, . . . , n — 1 at a cost /its*. The problem is to find a feasible flow that minimises 
the costs. 

Next, we will argue that we only need to look for optimal lot-sizing solutions 
that occupy a tree in the directed graph constructed above. We do not claim that all 
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Fig. 10.3. There exists an optimal solution that lives on a tree. 



optimal solutions have this structure, but there is always an optimal solution that has 
this structure, and that we only need to look for this solution. The special structure 
of this solution makes it possible to find it with dynamic programming. 

Proposition 10.1. ULS has an optimal solution (x,s) with the following special 
structure, illustrated in Figure 10.3, 

i) st~ix t = for t = l,...,n, that is, production occurs only when the stock is 
empty, 

ii) if x t > then x t — Y^i=t d i for some k > 0, that is, if production takes 
place during period t then the amount satisfies exactly the demands of periods 
t,...,t + k. 

Proof. Let (s, x, y) be an optimal solution occupying a minimial set of arcs. The 
result follows immediately from the no-wastage constraints 

s t _i + x t — d t + s t 

if we can show that the set of arcs on which (x, s) > is a tree in the flow-graph G 
constructed above. To prove this, it suffices to show that if 

x T , xe > 0, 

x T+ x, . . .,xe-i = (10.1) 
for t < 9, then sg-i = 0. (10.1) implies that 

s r = s T+ i + d T+ i > s T+ i > ■ ■ ■> S0-i. 
Let r := pg — (p T + h T + • ■ • + h$-i) and consider the feasible flow 







f- xe, 




= s t -f 


-x e it — t , . . . , 9 1), 


xe 


= 0, 




Xt 


= x t , 


(t^{r,6}), 


h 


= s t , 


{t£{r,...,9-l}, 


m 


= 0, 




m 


= Vt, 


(t ? 9). 
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Writing g(s, x, y) = X)"=i Pt x t + Y%=i + SLi 1 ^* s * for tne objective function, we 
have 

9(s,x,y) - g(s,x,y) = -(fe + rxg) > 0, 
where the last inequality follows from the optimality of (s,x,y). Therefore, 

r < -fo/x < 0. 

Next, consider the feasible flow 





= x T - se-i, 




St 


= S t -Sg-i (t = T,.. 


.,9-1), 


Xg 


= Xg + Sg-i, 




Xt 


= x t , {ti{T,9}), 




St 


= s t , (t<£{T,...,6- 


1}, 




_ JO if x T = s _i, 




Vr 


1 1 otherwise, 




Vt 


= Vt, (t^r). 





Then 



g(s,x,y) - g(s,x,y) 



rsg-i if x T > Sg-i, 

rSg-i-f T \ix T = Sg-i. 



Since we have already found r < 0, this contradicts the optimality and/or edge mini- 
mality of (s,x,y), because the flow (s,x,y) occupies one arc fewer than (s,x,y) and 
has a lower objective function. □ 

To simplify the algorithmic treatment of ULS, we eliminate the variables St from 
the problem. Using induction will now prove that the no-wastage constraints imply 
the relationship 

t 

i=l 

which allows to express the decision variables s t as a linear function of the decision 
vectors x,y. 

Lemma 10.2. Relation (10.2) holds for t = 1, . . . , n. 
Proof. The no-wastage constraint at time 1 

x 1 =d 1 +s 1 (10.3) 

yields si = x\ — d\, showing that (10.2) holds for t = 1. Assuming that (10.2) holds 
for t, the no- wastage constraint at time t + 1 

st + x t+ i = dt+i + s t+ i 
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implies 



t 

s t +i =s t + x t +i - dt+i = y^jxj - di) + (x t+ i - dt+i). 

i=l 

□ 



In eliminating the variables s, note that the constraint so = was already built 
into (10.3). The constraint s n = on the other hand is superfluous, since a production 
schedule with s n > can clearly not be optimal. 

Next, let us write the objective function in terms of the variables x,y. We have 

n n n— 1 t 

g(s,x,y) = ^2 Pt x t +^2f t y t + h t ^2(x t - d t ) 

1 = 1 t=l t=l 2 = 1 

n n n—1 

= ^2 c txt + ^2 ftyt - ^2 htdu =: ^ y )' 



t=i t=i t=i 



where we used the notation c t := pt + J27=t ^* and d it := J2j=i^j- Since the term 
— J2t=i htdit is just an additive constant, the minimiser of £(x, y) is the same as the 
minimiser of 



^2<hx t + ^2 ftyt- 
t=i t=i 



We now claim that the following problem is equivalent (in the sense that they 
have the same optimiser (x, y)) to the UCL problem, 

n n 

(P) H(n) = min V c t x t + V" f t y t 

x.y 

t=l t=l 

s.t. x t < My t (t = l,...,n), 

x e Rl, y £ B™, 

Vti < t 2 , x tl > = x tl+ i = ■■■ = x t2 -i < x t2 
=>x tl =d tl + --- + d t2 -i (10.4) 

We have already discussed why the omittance of the constraints so,s n = is 
okay, and why it is okay (without loss of generality) to impose the constraints (10.4). 
In addition, we have dropped the constraints s t > 0, since these are automatically 
satisfied when x satisfies (10.4). 

Note that (P) is no longer an IP formulation, but it has the advantage that lower- 
dimensional analogues can naturally be formulated by focussing on production and 
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consumption during periods 1, . . . , k only (in other words, by replacing n by k), 

k k 
(P fe ) H(k) = min ^ c t x t + £ f tVt 
x ' y t=i t=i 
s.t. x t < My t (t = l,...,k), 

xeR k + , y eB k 

Vii < t2, x tl > = = • • • = x t2 -i < x t . 2 

=> x tl = d tl H hdt 2 -i. 



Lemma 10.3. /iTey Observation] Let (x,y) be an optimal solution of (P) that 
satisfies the structure described in Proposition 10.1, and let t — max{# : yo — 1} be 
the last period during which production occurs in the schedule (x,y). Then 

i) x t = dtn, and 

ii) H(n) = H(t — 1) + f t + c t d tn (the "Bellman principle"). 

Proof. Part i) is an immediate consequence of the tree structure on which (s, x, y) 
lives. Furthermore, since St-i — 0, the production schedules periods 1, ... ,t — 1 is 
completely independent of the production schedule for periods t, . . . ,n. Therefore, 
(xi, . . . ,xt-i), (j/i, . . . ,y t -i)) must be a minimiser of (P t _i), and this shows ii). □ 

Algorithm 10.4. [Dynamic Programming for ULS] 

1. Initialisation: Set H(0) = 0, ro = n + 1, £ = 0, x, y = 0. 

2. Forward recursion: For k = l,...,n, compute 

H(k) = min (H(t - 1) + f t + c t d tk ) . 

t=l,...,k 

(Here we exploit Lemma 10.3.U).) 

3. Backward recursion: While 17 > ; repeat 

i) re+i = argmin^i,...^.! (H(t - 1) + f t + ctd^^^), 

Vr e+1 = 1, 
Hi) x Ti+1 d Ti+l T£ — i, 
iv) increment £ by 1. 



Theorem 10.5. Algorithm 10. 4 solves the ULS problem (P) in G(n 2 ) time. 

Proof. The correctness of Algorithm 10.4 is a direct consequence of Lemma 10.3. 
Since there are n iterations, each costing G(n) time because up to n values have to 
be compared with one another, the algorithm terminates in ff(n 2 ) time. □ 

Example 33. Consider the ULS problem with n = A, d = (2,4,5,1), p = 
(3, 3, 3, 3), h = (1, 2, 1) and f = (12, 20, 16, 8). 

• Using c t = p t + J27=t h i, we find c = (7, 6, 4, 3). 

• Using d it :=Ylj=i d j> we fi nd ( d n> d i2, d 13 , d 14 ) = (2, 6, 11, 12). 
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• H(l) = H(0) + h + Cl d n = 26. 

• Using the forward recursion, we find 



H{2) = min{# (0) + fi + Cl d 12 ,H(l) + f 2 + c 2 d 22 } 
= min{26 + c\d 2 , 26 + f 2 + c 2 d 2 } 
= min{54,70} = 54. 

• Continuing the forward recursion, we find 

H(3) = min{54 + Cl d 3 , 70 + c 2 d 3 , 54 + f 3 + c 3 d 3 } 
= min{89,100,90} = 89, 

H(4) = min{89 + Cid 4 , 100 + c 2 d 4 , 90 + c 3 d 4 , 89 + / 4 + c 4 d 4 } 
= min{96, 106, 94, 100} = 94. 

• Working backwards, we find H(A) = H(2) + f 3 + c 3 d 34 , so that y 3 = 1, 
x 3 — d 3 + d 4 = 6, and y 4 = 0, x 4 = 0. 

• Working backwards from t = 2, we find H{2) = H(0) + f\ + C\d\ 2 , so that 
yi = 1, x\ = d\ + d 2 = 6, and y 2 = 0, x 2 = 0. We have thus identified the 
optimal production plan x — (6, 0, 6, 0) with corresponding y = (1, 0, 1, 0). 



10.2. Binary Knapsack Problems. Recall from earlier parts of this course 
that the knapsack problem 

n 

max CjXj 

n 

s.t. ^ ajXj < b 
x e 7L\ 

can be relaxed by the integer knapsack problem with integer coefficients 

n 

(R) max ^ CjXj 

3 = 1 

n 

s.t. ^2[aj\xj < [b\ 

x e 7L\ 

In other words, the relaxed problem (R) is a knapsack problem in which the volumes 
aj and b are positive integers. The values Cj can be arbitrary real numbers. Integer 
knapsack problems can be solved via dynamic programming. Before we discuss the 
general case we consider the situation where x is a binary vector, that is, a single 
object is available of each type of item that can be packed into the knapsack. 
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Let thus cij 6 N™, fee N, and c £ K™ be given vectors. We shall see that the the 
0-1 knapsack problem 



(P) max cjXj 



3 = 1 

n 

s.t. ^ cijXj < b 

3 = 1 

x G B" 

can be solved in ff(nb) time. Note that it takes only ^(log&) bits to represent b in 
computer memory, so a running time of (?(nb) is actually exponentially large as a 
function of the input size. Nonetheless, for a moderate size of b, this approach is 
reasonable, as the complexity is only exponential in the size of b, and not in the size 
of a or cii. 

To set up a recursive hierarchy of binary knapsack problems, consider 

r 

(-Pr(A)) / r (A) = max cjXj 
i=i 

T 

s.t. cijXj < A 

3 = 1 

xeW. 

Then z = /„ (b) gives us the optimal value of our original knapsack problem. Further- 
more, the problems (P r (A)) are also binary knapsack problems but of smaller size. To 
arrive at a recursive relationship between the values / r (A), we distinguish two cases: 

• If x* = 1, then the choice of the remaining variables x±, . . . ,x r _i must be 
optimal for (P r _i(A — a,.)), and hence, 

/ r (A) = c r + f r -i(X -a r ). 

• If x* = 0, then the choice of the remaining variables x\, . . . ,x r _i must be 
optimal for the problem (P,._i(A)). In other words, we then have 

frW = /r-l(A). 

How do we know in which case we are? Just compare the two function values we 
would get for / r (A) and pick the one that produces the larger value! Thus, we arrive 
at the recursion 

/ r (A) = max(/,._i(A),c r + / r _i(A - a,.)). 

To initialise the recursion, we set the obvious boundary values 

/ r (0) = 0, (r = l,...,n), 
/o(A)=0, (A = 0,...,6). 



Algorithm 10.6. [Binary knapsack with integer coefficients] 
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Initialisation: Set / r (0) = 0, (r = 1, . . . , n), /o(A) =0, (A = 0, . . . , b). 
Forward Recursion: 

For r = 1, . . . , n, repeat 

/r(A) = / r -i(A), (A = 0,...,a r -1), 

/ r (A) = max(/ r _i(A),c r + / r _i(A - a r )), (A = a r , . . . ,6). 

end. 

Backward Recursion: Set X = b, r = n, x = 0. 
WMe r, A > 0, repeat 

if /r(A) = c r + / r _i(A - a r ) 
;c r = 1 ; A = A — fj r 

end 

r = r — 1 

end. 



There are &(nb) inner loops both for the forward and backward recursion, each 
costing 0(X) time. Therefore, the complexity of this algorithm is ff(nb). The com- 
putations of the backtracking procedure can be avoided altogether by keeping track 
of the indicator variables 



p r (X) = 



if /r(A) = / r _i(A), 

1 if / r (A) = c„ + / r _i(A - a n ) 



during the forward iteration. This comes at the obvious expense of increasing the 
memory needs because the values of p r (X) need to be stored. 

Example 34. Let us apply Algorithm 10.6 to the 0-1 knapsack problem 



max 10a; i + 7x 2 + 25x 3 + 24x 4 
s.t. 2x\ + x 2 + 6x 3 + 5x 4 < 7 
x G B 4 . 



• Starting the forward recursion, /i(l) = /o(l) = 0, Pi(l) = 0, 

• /i(2)=max(/o(2),ci+/ (2-ai))=max(0,10 + 0) = 10, pi(2) = I, 

• and likewise, /i(3) = • • • = /i(7) = 10, pi(3) = • • • = pi(7) = 1. 

• Next, / 2 (1) =max(/ 1 (l),c 2 + /i(l - o 2 )) - max(0, 7 + 0) = 7, p 2 (l) = 1, 

• /2(2)=max(/ 1 (2),C2 + /i(2-a 2 ))=max(10,7 + 0) = 10, p 2 (2) = 0, 

• / 2 (3) = max(/ 1 (3),c 2 + A(3 - a 2 )) = max(10,7+ 10) = 17, p 2 (3) = 1, 

• and likewise, / 2 (4) = • • • = / 2 (7) = 17, p 2 {A) = ■■■ = p 2 (7) = 1, 

• Continuing this way, we find the values 
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Back-tracking is now easy, as we have computed the indicators p r (A) : 

• Pi(7) = 1, thus x\ = 1. 

• Therefore, we next have to check p?,(7 — 04) = pa(2) = 0, so Xg = 0, 

• so we check ^2(2) = 0, showing that x\ = 0, 

• and finally, pi(2) = 1, so that x\ = 1. 

Thus, x* — (1, 0, 0, 1) is an optimal solution. We check that x* is indeed feasible, 

a T x* = 2 + 5 < 7, 
and achieves the promised objective value, 

c T x* = 10 + 24 = 34. 



10.3. Integer Knapsack Problems. Next, we extend our method to general 
integer knapsack problems in which multiple copies of the same item are available for 
packing, 

n 

max CjXj 

n 

s.t. ajXj < b 

x e U\. 

The coefficients 6, a 3 £ N (j = 1, . . . , n) are again assumed to be positive integers. 

We can again vary the right-hand side A = 0, . . . , b and the number of variables 
r = 0, . . . , n by considering all of the following integer knapsack problems, 

r 

g r (X) := max ^ CjXj 

r 

S.t. OjXj < A 

x e 7T + . 
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Clearly, our original knapsack problem is (P„(6)), so that its optimal value is given 
by g n {b). Furthermore, the following boundary values are obvious, 

g r (0) = 0, (r = 0,...,n) 
S0(A) = O, (A = 0,...,6). 

Note that at most copies of item r fit into a knapsack of volume A. Using 
the Bellman principle, and distinguishing the cases 



x* = 0, . . . 

we thus find the recursion formula 

<7 r (A) = max |te r + fif r -i(A — ta r ) : t = 0, . . . , 

However, each of these iterations takes 0(\b/aj\) to compute, leading to an algorithm 
with an overall complexity of 



^n&max [b/a 3 ^j = (nb 2 ) . 



We can do better than this: Observe that if x* = 0, then the vector {x\, . . . , x*_ 1 ) 
must be optimal for the problem (P r _i(A)), so that g r (X) = g r -i(X). If on the other 
hand x* > 1, then the vector (x\, . . . , x*_ x ,x* — 1) must be optimal for (P r (A — a r )), 
so that <7r(A) = c r + g r (X — a r ). Therefore, we can use the recursion 

g r {\) = max(g r _i(A),c r + .g r (A - o r )), 

which takes time to compute. We need to compute &(nb) such recursive steps, 
so we obtain a 6\nb)-\\me algorithm: 

Algorithm 10.7. [Integer knapsack with integer coefficients] 

Initialisation: Set g r {0) = 0, (r = 1, . . . , n); 50(A) =0, (A = 0, . . . , b). 

Forward recursion: 

for r = 1, . . . , n, repeat 

g r (X) = g r -i(\), A = 0, . . . , a r — 1 

g r (X) = max(g r _i(A),c r + g r (X - a r )), (A = a r , . . . ,b) 

end. 

Backward Recursion: Set X = b, r = n, x = 0; 
while A, r > 0, repeat 
if 9r{X) = SV-i(A) 

r = r — 1, 
elseif g r (X) = c r + g r (X — a r ) 

X vf* — X y 1 . 

A = A a r , 

end. 

For the purposes of the backward sweep, it is again convenient to collect the 
values of the following variables, 



Pr{X) 



1 if g r (X) = c r +g r (X - a r ), 
otherwise. 
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To back-track, we have to start checking p n (b): If p n (b) = 1 then x* > 1 and we must 
check p n (b — a n ) to see whether x* n > 2 etc.. If on the other hand p n (b) — 0, then 
x* = and we must check p„_i(6) to see whether x* n _ x > 1 etc.. 

Example 35. Consider the integer knapsack problem 

max7xi + 9x2 + 2x3 + 15x4 
s.t. 3xi + 4x 2 + x 3 + 7x 4 < 10 
x G 1\. 

Applying Algorithm 2, we find the following table of values for g r (X), p r (ty, 
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Back-tracking: 

• p 4 (10) = 0, so x* = 0. 

• p 3 (10) = 0, so xl = 0. 

• p 2 (10) = 1, so x* 2 > 1. 

• P2(10 — a 2 ) = p 2 (6) = 0, so x 2 ^ 2, and hence, x 2 = 1. 

• pi(6) = 1, so xl > 1. 

• pi(6 - ai) = pi(3) = I, so x* > 2. 

• pi(3 — ai) = pi(0) = 0, so x\ ^ 3, and hence, x* = 2. 

We have found that x* — (2, 1, 0, 0) is an optimal solution. 

Knapsack problems with integer coefficients can also be reformulated as a longest 
(or shortest) path problem. Consider again the problem 

n 

max CjXj 

3=1 

n 

s.t. ajXj < b 

3=1 

x G Z™ 

with b, dj G N (J = 1, . . . , n), and construct an acyclic digraph D = (V, A) with 

• node set V = {0, 1, . . . , b}, 
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Fig. 10.4. The knapsack problem from Example 35 drawn as a directed graph. 



• arcs set 

A = {(A, A + a,j) : A e Z + s.t. A < b — a,-, j = 1, . . . ,n} 
U{(A,A + 1) : A<6-1}, 

• arc weights Cj for the arcs (A, A + a,j), 

• and arc weights for the arcs (A, A + 1). 

Then <m(A) is the length of a longest path from node to A. 

For example, the diagraph corresponding to the knapsack problem from Example 
35 is as drawn in Figure 10.4. We already saw that the longest (or shortest) path 
problem can be solved in polynomial time, using network flows. Note however that 
the digraph corresponding to a knapsack problem contains exponentially many arcs 
as a function of the input size. 

11. Lagrangian Relaxation. Consider the integer programming problem 

(IP) z = ma,xc T x 

s.t. Ax < a 
Dx <d 

x e z%, 

where Ax < a is a benign set of constraints (e.g., totally unimodular) in the sense 
that 

T 

maxc x 
s.t. Ax < a 

x e Z™ 

would be easy to solve, and where Dx < d is a set of m malicious constraints that 
render (IP) intractable (e.g., connectivity constraints in the TSP). 

For such problems we will now derive a family of relaxations that can generate 
stronger bounds than LP relaxations. Consequently, branch- and-bound systems built 
around these relaxations are usually more efficient than an LP based approach. We 
start by writing (IP) in the slightly more general form 

(IP) z = maxc T x 

s.t. Dx@d, 

x e x, 
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where @ can stand either for "<" or "=", and where X is a feasible set of "benign" 
type. Let us write J for the set of indices that correspond to inequality constraints 
among the system Dx@d, and £ for the set of indices that correspond to equality 
constraints. 

Definition 11.1. A Lagrangian relaxation of (IP) is a problem of the form 
(IP{u) ) z(u) = max{c T x + u T (d - Dx) : x e X} 
where u 6 R m is a fixed vector Lagrange multipliers, and where Ui > for ie/. 

Proposition 11.2. (IP{u)) is a relaxation of (IP). 

Proof. Firstly, the feasible region of (EP(tt)) contains that of (IP), 

X D {x G X : Dx@d}. 

Secondly, for all x feasible for (IP) the objective function of (EP(u)) is at least as large 
as that of (IP), 

c T x + u T (d — Dx) = c T x + ^ Ui(di — Di t -x) > 0, 

where we have used that for feasible x, di — D^-.x > for alH £ J 1 and di — D it -x = 
for all ieS.U 

Obviously, there exist infinitely many Lagrangian relaxations (IP(u)) to choose 
from. So how should we choose the vector ul Since (IP(u)) is a relaxation of (IP), 
its optimal solution yields an upper bound z(u) on the optimal objective value z of 
(IP). It is therefore natural to seek the least upper bound of this kind. This leads to 
the Lagrangian dual problem 

w LD = min{z(w) : u £ K m , u t > OVi G J}. 

As with any relaxation technique, Lagrangian relaxation can sometimes yield an op- 
timal solution to the original problem (IP): 

Proposition 11.3. Ifx(u) is optimal for (IP(u)) and satisfies Dx(u)@d and 

{Dx(u))i=di Mi e J s.t. Ui> 0, (11.1) 

then x(u) is optimal for (IP). 
Proof. 

Propll.2 (n -|\ 

z < z(u) ^ c 1 x(u) + u 1 (d- Dx(u)) = c L x(u)<z. 

□ 
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11.1. Lagrangian Relaxation of the UFL. Lagrangian relaxation is a partic- 
ularly powerful approach in the context of the uncapacitated facility location problem 
we encountered in Section 1, 

(IP) z = max Ci i Xi i ~ Y UVi 

ieM jeN jeN 

S.t. Y X ij = 1 (* G M ) 

jeN 

x tJ -yj<0 (ieMJeN) 




where M is the set of customer locations, N is the set of potential facility locations, fj 
are the fixed costs for opening facility j, and where we replaced the original servicing 
costs Cij with —Cij to turn the problem into a maximisation problem. 

"Dualising" the demand constraints ^2j eN x^j = 1, we find the Lagrangian relax- 
ation 

(IP(u)) z(u) = max Y ( Ci i ~ u ^ ~ fiVi + u i 
ieM jeN jeN ieM 

s.t. Xij - yj < (i e M, j G N) 




Now note that the constraint that linked the different facility locations to one another 
has been subsumed in the objective function, so that (IP(w)) decouples, 

jeN ieM 
where Zj (u) is the optimal solution of the following problem, 

(IPj (u)) Zj (u) = max ^ ( c y ~ u i) x H - fiVi 

ieM 

s.t. Xij - yj < (i e M) 
x^ >0(j£ M), yj e B. 



Furthermore, (IPj(u)) is easily solved by inspection: 

• If yj = 0, then Xij — for all i, and the objective value is 0. 

• If yj = 1, then all clients i for which Cjj — u» > will be served, and the 
objective value is ^2 i( z M max(0, Cjj — Uj). 



Therefore, Zj(ti) = max^O, X)ieM max (0, — Ui) — f^j . 



11.2. Lagrangian Relaxation of the STSP. Let us now derive a Lagrangian 
relaxation of the symmetric TSP problem. We saw that this problem can be modelled 
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as follows, 



(IP) z = min E c e x e 

S.t. E X e = 2 ( i( ^ V ) 

ee8(i) 

£ x e <\S\-l (S C V s.t. 2 < |5| < |V| - 1) 

e£E(S) 

where V is the set of nodes, E is the set of (undirected) edges, c e is the cost of edge 
e, (5(i) C _B is the set of edges incident to node i, and E(S) is the set of edges in E 
that have both end points in 5. 



Lemma 11.4. Half the subtour elimination constraints 

£ *e<|S|-l 

eeE(S) 

are redundant. 

Proof. For any x feasible for the LP relaxation of (IP) we have 



\s\- e ^^EE 1 ^ £ 

ies ee 



2 ^ ■ , £<= 

eEE(S) lES eeS(i) e£E(S) 



2 

ee<5(S,S<=) 



where S(S, S c ) is the set of edges in E with one end point in S and the other in 
S c := V \ S. Since 5(S, S c ) = 8(S C , S), we now have 



£ x * = \ E * e = is c i- E 

eeE(S) eeS(S,S") eeE(5 

and hence, J2 e eE(s) < l^l - 1 ^> J2 e eE(S") x e ^ \ S °\ ~ L n 



Eliminating all subtour elimination constraints for which 1 £ S, dualising all 
degree constraints 

£ *e = 2 

e€«(») 

except for i = 1, and introducing the (redundant) constraint 

eS-E 

obtained by summing all degree constraints, we arrive at the following Lagrangian 
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relaxation of the STSP, 

(IP(u)) z(u) — min (c e — u; t — Uj)x e + 2 Ui 

e=(ij)£E ieV 
X e = 2 

e€<5(l) 

\ s \ _ x ' V5 c v s - t - 2 < |S| < |V| - 1. 1 ^ 5 

eGE(S) 

= n 

eS-E 



Definition 11.5. 4 1-tree m an undirected graph G — {V,E) is a subgraph 
that consists of the union of two edges adjacent to node 1 and a spanning tree on the 
remaining nodes. 

Lemma 11.6. x s Bl B l is feasible for (IP(u)) if and only if the subgraph G x := 
(V, {ee E : x e = 1}) is a 1-tree in G = (V,E). 

Proof. The constraint J2ee6(i) Xe = 2 guarantees that exactly two edges are 
incident to node 1 in the subgraph G x := (V, {e € E : x e = 1}. The constraints 
See-E(S) Xe — 1^1 — 1 guarantee that when node 1 is removed, then there is no cycle 
left in the subgraph G x . The constraint J2eeE x e = n guarantees that there are 

n-2=\V\{l}\-l (11.2) 

edges in G x \ (5(1), and since G x \ 6(1) does not contain any cycles, (11.2) shows that 
it must be a tree, and hence, G x is a 1-tree. On the other hand, it is straightforward 
to check that the incidence vector a; of a 1-tree satisfies the constraints. □ 



Combining the insight of Lemma 11.6 with the results of Section 8, it follows that 
(IP (it)) can be solved via a greedy approach: Build a maximum weight spanning tree 
on the nodes {2, . . . , n} and add the two heaviest edges incident to node 1. 



Example 36. Let us now look at a numerical example and consider the STSP 
on 5 nodes with edge cost matrix 

~- 30 26 50 40 

30 - 24 40 50 

26 24 - 24 26 

50 40 24 - 30 

40 50 26 30 - 

Note that the Lagrange multipliers u are unrestricted, as we have dualised equality 
constraints. Therefore, 



u=[0 -15 0] 
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is a legitimate choice. Writing Cij := Cij — Ui — Uj, we obtain the revised edge cost 
matrix 





30 


41 


50 


40 


30 




39 


40 


50 


41 


39 




39 


41 


50 


40 


39 




30 


40 


50 


41 


30 





To build a minimum weight 1-tree, we 

i) find a minimum weight spanning tree on the nodes {2,3,4,5} 
ii) and add two edges of minimum weight incident to 1. 

In an earlier sections we have seen that the greedy algorithm solves part i): Order the 
edges according to increasing edge weight and add them in this order as long as they 
don't create a cycle. Our edges are ordered as follows, 

(4,5),(2,3),(3,4),(2,4),(3,5),(2,5). 

We choose (4,5), (2,3), (3,4) in this order and then stop because adding any of the 
other edges would create a cycle. Further, part ii) is trivial to solve, as (1, 2) and (1, 5) 
are the cheapest edges incident to node 1. Therefore, {(1, 2), (1, 5), (4, 5), (2, 3), (3, 4)} 
is an optimal 1-tree. The corresponding incidence vector is actually is STSP-feasible 
and corresponds to a Hamiltonian tour 

l->2->3-»4->5->l. 

By Proposition 11.3, it follows that this is an optimal STSP-tour. 

11.3. Solving the Lagrangian Dual. In order to make Lagrangian Duality a 
useful tool, we need an algorithm that can solve the Lagrangian Dual. Let us continue 
to consider IPs of the form 

(IP) z = maxc T x 

s.t. Dx < d 

x e X := {y e Z™ : Ay < b}, 

where Dx < d are m constraints that make the problem difficult, while optimising 
over X is easy. In Section 11 we also considered equality constraints Dx = d, but 
we will not treat them separately anymore, as they can easily be subsumed in the 
standard formulation (IP) in the form 

Dx < d, 
-Dx < -d, 

or else they can be treated separately by dropping the bound constraints on the 
associated Lagrange multipliers. 

The Lagrangian dual of (IP) is 

(LD) wld = minz(u) 
s.t. u 6 M+ , 
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where 



(IP(u)) z(u) — maxc T a; + u T (d — Dx) 
s.t. x e X. 

We saw that (LD) is a dual of (IP), so that z < wld- Solving (LD) is thus a good way 
of generating a hopefully quite tight bound on z, but note that u i->- z(u) is generally 
a nonlinear function, and we do not know yet how to solve such problems. 

To get around this problem we will have to understand the structural properties 
of the map u i-> z(u). This structure is best revealed by asking the question as to 
whether Lagrangian duality yields stronger bounds than LP relaxation. To answer 
this question, let us first assume that X contains a large but finite number of points 
X = {jW,...,^ 1 }. Then 

wld = niin z(u) 

u>0 

= min{max{c T x [tl + u T (d - Dx [t] ) : t = l,...,T}\ 

= min r\ 

(u,t))ei m + 1 

s.t. T]>c T xW+u T (d-DxW), (t=l,...,T) 
u > 0. 

The latter problem is an LP instance with decision variables {u,rj). Taking the LP 
dual, we obtain 



wld = max 



T 

E 



s.t. j2^( Dxlt] - d ) < ° 



t=i 

T 

t=i 

A* > 0. 



Since 




A*>0, 5Z/*t = l| = conv ({xW, . . . , x^}\ 
t=i J 



is the convex hull of X, we thus arrive at the following result: 
Theorem 11.7. 



wld = maxc T x 



s.t. Dx < d 

x G conv(X). 
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We remark that this result still holds true for arbitrary feasible sets of the form 



X = {xeZ r l: Ax< 6}, 

not just in the finite case. 

Corollary 11.8. 
i) If {x G R™ : Ax < b} is an ideal formulation of X , then 

conv(X) = {i€l": Ax < b, x > 0}, 

and hence, wld coincides with the bound produced by the LP relaxation of 
(IP), 

wld = max{c T a; : Dx < d, Ax < b, x > 0}. 

ii) If {x G K" : Ax < b} is not an ideal formulation of X then 

conv(X) c{i£ R" : Ax < b] 

and wld is generally a tighter bound than the one given by the LP relaxation 
of (IP)- 

Let us remark that in situation ii) one is rewarded with a stronger bound, but this 
comes at a price since the fact that {x G R™ : Ax < b} is not an ideal formulation of 
X means that (IP(m)) may not be that easy to solve. Situation i) is not rare, because 
max{c T x : x G X} is an easy problem mainly when {x G M™ : Ax < b} is an ideal 
formulation of X. Solving the Lagrangian dual is nonetheless interesting in this case, 
as solving the LP relaxation directly may be too costly. For example, in Section 11 
we found that a Lagrangian relaxation of the STSP is given by 

(IP(it)) z(u) = min (c e — Uj — Uj)x e + 2 Uj 

e=(ij)£E ieV 

s.t. x e = 2 

e€«(l) 

^2 x e < \S\ - 1, VScF s.t. 2<\S\<\V\-1, l^S 

eeE(S) 

^x e = n 

.tgb! £ I. 

Solving the LP relaxation directly would be difficult, as it contains an exponential 
number of constraints. However, solving (IP(u)) via the greedy algorithm is simple, 
and we will next see that when (IP(u)) can be solved efficiently, then the Lagrangian 
Dual can in general also be solved reasonably efficiently using the subgradient algo- 
rithm. 

The above discussion revealed that z(u) is a piecewise linear function 

z(u) = t _max r (c I V*] +u T {d~Dx^)). (11.3) 
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Furthermore, this function is convex, since the linear functions u \— > c T x^ + u T (d — 
Dx 1 ^) are convex and the pointwise maximum of a set of convex functions is again 
convex. Note however that z(u) is not diffcrentiable at the "breakpoints" where the 
maximum in (11.3) is achieved by more than one index, and that to solve (LD) we 
need to minimise z(u) over u > 0. Convex nondifferentiable functions can often be 
reasonably well solved by the subgradient algorithm. Before we can explain the algo- 
rithm, we need to understand the notion of a subgradient. 

Lemma 11.9. Let f : W n — > R be a convex function with gradient 7 = V/(u) at 
u E M. m . Then the first order Taylor approximation satisfies 

fi^+^iv-u) <f(v) 

for all v eR m . 

Proof. By definition, / is convex if for all u, v £ K m and A € [0, 1], 
/ (Xv + (1 - A)u) < Xf(v) + (1 - A)/(u). 

Therefore, 

m + fiu+x{v - u)) - m <m, 

and taking the limit A — > yields the result. □ 

The notion of subgradient now generalises the property we just observed at points 
where / is nondifferentiable: 

Definition 11.10. Let u e W n and let f : M m ->■ M be a convex function. A 
subgradient of f at u e R m is a vector 7 e R m such that 

/( U )+7 T («- U ) < f(v) 

for allv E R m . 

We remark that if / is diffcrentiable at u, then V/(u) the unique subgradient of 
/ at u. Further, any convex combination of subgradients is a subgradient. The set of 
subgradients of / at u is thus a convex set V/(u) called the subdifferential of / at u. 
To apply the subgradient algorithm in solving the Lagrangian dual, we will need to 
be able to compute elements of Vz(u). The next result yields a mechanism to do that. 

Lemma 11.11. Let u e K m and 

t* = arg max jcV* 1 + u T {d - Dx [t] ) : t = 1,...,t}. 

Then d - Dx^ is a subgradient of z(u) at u. 
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Proof. For all ueT, 

z{v) = t max^cV^ + v T (d - Dx&)) 
>c T x^+v T (d-Dx^) 

= c T x^ + u T (d - DxW) + (v- u) T (d - Dx^) 
= z(u) + (v-u) T (d- D/l). 

By now we have all the tools in place to discuss the subgradicnt algorithm: 

Algorithm 11.12. [Subgradient Algorithm for Solving (LD)] 

1. Initialisation: Choose initial Lagrange multipliers u'°' > via educated guess. 

2. For k = 0,1,2, ... repeat, 

i) find a maximiser x(u^) of the Lagrangian relaxation (IP(u^)) corre- 
sponding to the Lagrange multipliers u^, 
ii) determine a step length multiplier \i k > 0, 
Hi) for i = 1, . . . , m, set 

uf +l] = max( U f ] - n k {d - Dx(u^ k %, 0) 

end. 



The motivation of the algorithm is very simple: in each iteration of the main loop 
the Lagrange multiplier vector is improved by correcting it in a direction that makes 
the objective function z(u) decrease. Note the built-in safeguard mechanism that pre- 
vents individual components of vJ k+1 ^ to become negative. For Lagrange multipliers 
corresponding to equality constraints Dx = d, this safeguard is of course unnecessary. 
The implementation of the subgradient algorithm is very simple too, but the difficulty 
lies in choosing the step length multipliers fik, as the following theorem shows. 

Theorem 11.13. 

i) If J2 k [iu — > oo and \ik — > as k — »• oo, then z{uV^) — > wld- 
H) If Mfe = PoP k for some fixed p G (0, 1), then z(W fe l) — > wld if Po an d P are 

sufficiently large. 
Hi) if w > wld and 

e k x (z(u[ fc l) - w) 
^ k = \\d - Dx(uW)\\ 2 ' 

where e k G (0,2) for all k, then either z(u[ fe l) — > wld for k — > oo, or else 
w > z(u^ ) > wld occurs for some finite k. 
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Example 37. Let us continue the the discussion of Example 36. The dualised 
constraints were the degree constraints 

*e = 2 (i e V). 

eeS(i) 

Since these are equality constraints, the Lagrange multipliers u are unconstrained, 
and the updating rule is 

tt [^=«{*]+ Mfc (2- ]T *„(«[*])). 

e£<5(») 

TTie STSP being a minimisation problem rather than a maximisation problem, we have 
to replace all mins by maxes, lower bounds by upper bounds and so forth. We will use 
rule Hi) of Theorem 11.13 to update Lagrange multipliers, using = 1 for all k. The 
step length multiplier is therefore given by 

w — z(u^) 

where w is a lower bound on wld, and where we had to invert the sign of fik because 
(LD) is a maximisation problem. 



Initialisation: We apply the greedy heuristic and find the tour 1^2^3^4— > 
5 ->• 1 of length 14.8. 

As no lower bound w on wld is known at present, we use the primal (upper) 
bound w = 148 instead and see how far this allows us to decrease z(u^), knowing 
that at a later time we will probably have to replace w by a smaller value, as 148 is 
usually not a lower bound. 

Iteration 0: Starting with = [0,0,0,0,0], the revised costs are given by 

-[0] _ _ [0] _ [0] _ 

C ij — C ij U i U j ~ C ij- 

Solving the associated min-cost 1-tree problem, we find the optimal 1-tree with 
edge incidence matrix 

-110 0" 

- - 1 

- - - 1 1 , 

- - - - 



leading to the objective value z(u^) = 130. 
The subgradient at x(u^) is 

[(2- *e(« M ))<=i,...,J = [0,0,-2,1,1], 
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[Sy (UN)] = 



and as /zo = (148 — 130)/6 = 3, we find 

U W = u [o\ + 3 . [0,0,-2,1,1] = [0,0,-6,3,3]. 



Iteration 1: The new cost matrix is 



Jill 



- 30 32 47 37 

- - 30 37 47 

- - - 27 29 

- - - - 24 



The optimal 1-tree is found as 



[ Xij {uM)] 



110 
-10 

--10 

- - - 1 



leading to the objective value z(u^) = 143 + 2 J2i u ?' = 143. 
Updating the Lagrange multipliers, we obtain 



„m = „[.] + M£_iH . [0,0,-1,0,1] 



Iteration 2: The new cost matrix is 



- 30 34.5 47 34.5 

- - 32.5 37 44.5 

- - - 29.5 29 

- 21.5 



The optimal 1-tree is found as 



M« [2] )] = 



- i 



1 

1 

- 1 

- - 1 



leading to the objective value z(u^) = 147.5. 

This means that 147.5 is a lower bound on the optimal value z of the STSP. But 
since [cij] is integer valued, this implies 

148 = [147.5] < z < w = 148, 

which shows that the greedy tour 1^2^3— ^4^5^1 was STSP-optimal! 
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11.4. Choosing a Lagrangian Dual. Many problems have several reasonable 
Lagrangian Duals. In this case it is worthwhile thinking about the advantages and 
disadvantages of the different formulations before starting any calculations. 



Example 38. Consider the generalised assignment problem 

n m 

(IP) z = max ^2 

CijXij 

j=i i=i 

n 

s.t. y^Zjj < 1 (i = l,...,m), 

m 

^dijXij < bj (j = l,...,n) 

i=l 

x G B mx ". 

In this case we have multiple choices of a Lagrangian dual: 
1. Dualising both sets of constraints 



(IP{u)) z(u) = max ^ X^ c ^ _ Ui ~ a ii v i) x io + X! Ui + X! 

J — 1 2—1 i=l J = l 



m n 



s.t. x G 



If is certainly easy to solve this IP, as the remaining feasible set X — B mx " has the 
ideal formulation {x G R m : < Xi < 1 Vi}. #?/ Theorem 2, solving the Lagrangian 
dual min„>o z{u) yields the same result as the LP relaxation of (IP). However, it might 
still be better to solve the Lagrangian dual, as (IP{u)) can be solved by inspection: 



x ij ~ 



1 ifcij - Ui - aijVj > 
otherwise. 



2. Dualising only the first set of constraints 

n m m 

(IP{u)) z(u) = max ^ ^(cjj - Ui)Xij + u { 

j — l i—l i—1 

m 

s.t. ^2 ' U = 1, ■•-,") 



ifere the remaining feasible set 



X = jz G B mx " : ^ a^xy < 6 j; j = 1, . . . , raj 



is not of the "easy" type, as 



r(X) C |i£ K™ : < ^ < 1, i = 1, . . . , m, ^ ( 



com | X ) c < .i fc A" : U .r, - 1 m, a^Xy <bjj = l,...,n 
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is generally a strict inclusion. Consequently, Wld may be a strictly tighter bound than 
the bound obtained from the LP relaxation, but the Lagrangian dual is more difficult 
to solve. 

3. Dualising only the second set of constraints The advantages and disadvantages 
of this relaxation are similar to the first case. 
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